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ABSTRACT 


System  identification  concerns  the  mathematical  modeling  of  a  system  based  upon 
its  input  and  output.  It  allows  the  development  of  a  mathematical  description  when  all 
that  is  available  is  the  result  of  a  process  or  the  output  of  a  system  and  not  the  process 
or  system  itself. 

The  purpose  of  this  thesis  is  to  develop  algorithms  for  modeling  systems  as 
autoregressive-moving-average  processes  using  the  method  of  instrumental  variables,  a 
modification  of  the  ordinary  least-squares  technique,  and  a  multichannel  method  based 
upon  processing  the  input  and  output  data  by  separate  infinite-impulse-response  filters. 
The  methods  developed  are  tested  by  computer  simulation  using  several  second  and 
third-order  test  cases  and  the  results  are  presented. 
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I.  INTRODUCTION 


A.  SYSTEM  IDENTIFICATION  BASICS 

System  identification  concerns  the  modeling  of  systems  as  sets  of  mathematical 
equations  based  upon  the  input  and  output  of  the  system.  [Ref.  1:  pp.  3-6).  It  allows 
a  model  to  be  developed  when  all  that  is  available  is  the  result  of  a  process  or  the  output 
of  a  system  and  not  the  process  or  the  system  itself.  System  identification  is  an  impor¬ 
tant  area  of  study.  Solution  of  the  modeling  problem  offers  many  alternatives  for  the 
continued  study  of  the  system.  Among  these  are: 

•  Nondestructive  analysis  of  the  system. 

•  Simulation  studies  using  the  model. 

•  Easy  adaptation  of  the  model  to  changing  system  environment. 

•  Spectral  analysis  of  the  system. 

Modeling  can  simulate  the  system's  operation  at  a  fraction  of  the  cost  of  actual 
system  operation.  Complex  operations  not  possible  with  the  actual  system  for  fear  of 
damaging  it  or  personal  injurs  can  be  simulated.  This  can  expose  howr  the  system  will 
operate  in  adverse  conditions  not  normally  experienced.  In  speech  processing,  modeling 
the  speech  process  has  the  potential  for  significantly  reducing  the  amount  of  information 
necessary7  to  store  in  order  to  reproduce  the  speech. 

The  modeling  process  shown  in  Figure  1  on  page  2  assumes  the  unknown  system's 
input  and  the  output  data  are  available  for  processing.  In  many  cases,  if  the  system's 
input  is  unknown  or  data  is  not  available,  a  white  noise  input  can  be  used  in  its  place. 
The  modeling  process  uses  the  input  and  output  data  to  find  a  set  of  parameters  which 
closely  approximate  the  operation  of  the  system.  The  better  the  identification  technique, 
the  more  closely  the  model  follows  the  performance  of  the  actual  system. 

Many  types  of  models  are  available.  This  thesis  investigates  a  linear  parametric 
model  that  can  be  described  by  difference  equations.  This  type  of  model  lends  itself  well 
to  simulation  on  a  digital  computer.  The  frequency  characteristics  of  the  system  deter¬ 
mined  from  the  parameters  of  these  types  of  models  are  more  accurate  than  wrhat  can 
be  determined  from  classical  means  such  as  FFTs.  This  is  because  classical  methods  use 
windows  w'hich  assume  data  beyond  their  extent  is  zero  [Ref.  2:  p.  173J.  This  is  not  a 
realistic  assumption.  Models  in  this  category  include  the  moving-average  (MA)  model, 
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Figure  i.  System  identification  problem 

the  autoregressive  (AR)  model,  and  the  autoregressive-moving-average  (ARMA)  model. 
In  the  frequency  domain,  MA  processes  are  characterized  by  sharp  nulls  and  smooth 
peaks  and  AR  processes  are  characterized  by  smooth  nulls  and  sharp  peaks.  ARMA 
processes  have  sharp  peaks  and  sharp  nulls  [Ref.  2:  p.  173).  An  advantage  of  the  MA 
process  is  its  inherent  stability.  An  advantage  of  the  AR  process  is  the  large  number  of 
algorithms  already  available  for  modeling  systems.  An  advantage  of  the  ARMA  process 
is  that  it  uses  far  fewer  parameters  than  either  the  MA  or  AR  process  alone  to  model  a 
system.  This  satisfies  the  general  requirement  to  reduce  the  complexity  of  the  model. 

In  addition  to  a  large  variety  of  models,  there  are  two  processing  modes:  block  and 
sequential. 

Block  processing  uses  a  fixed  length  block  of  data  in  the  parameter  estimation 
process.  It  ignores  data  before  and  after  the  block.  This  is  not  a  real  time  processing 
method  because  all  data  must  be  available  before  processing  can  start.  Block  processing 
generally  involves  inversions  of  data  matrices  whose  sizes  are  on  the  order  of 
(.V  +  .\f)  x  (.V  +  .\f)  where  .V  is  the  order  of  the  AR  process  and  M  is  the  order  of  the 
MA  process. 


Sequential  processing  uses  new  data  to  update  the  parameter  estimations.  It  starts 
by  initializing  an  estimate  of  the  inverse  of  the  data  covariance  as  as  a  diagonal  matrix. 
It  uses  each  new  data  point  to  update  this  matrix.  Then  it  updates  the  parameter  esti¬ 
mates  using  the  updated  inverse  data  covariance  matrix.  It  is  a  real  time  method.  The 
algorithm  to  implement  the  sequential  processing  method  is  generally  more  complex 
than  the  block  method  but  less  computationally  intensive  because  the  matrix  inversions 
are  not  required. 

This  thesis  concerns  only  systems  represented  by  discrete  time  data  uniformly  sam¬ 
pled  at  a  sufficient  rate  to  meet  the  Nyquist  criteria. 

The  work  in  this  thesis  assumes  that  the  input  data  is  a  wide-sense  stationary  ran¬ 
dom  sequence.  Tests  of  the  algorithms  used  a  pseudorandom  Gaussian  input  with  a 
mean  of  zero  and  a  variance  of  one. 

B.  PROBLEM  STATEMENT 

The  purpose  of  this  thesis  is  to  develop  algorithms  for  modeling  systems  as  ARMA 
processes  using  the  method  of  instrumental  variables  (IV)  and  a  multichannel  approach. 
Tests  of  the  methods  will  be  conducted  to  determine  the  accuracy  of  their  results  and  the 
speed  with  which  they  converge. 

The  IV  approach  is  a  modification  of  the  method  of  ordinary  least  squares.  This 
approach  is  developed  first  as  a  block  processing  case  and  then  converted  to  a  sequential 
processing  case.  Tests  are  conducted  of  only  the  sequential  processing  case. 

Using  a  multichannel  scheme  allows  the  input  and  output  data  of  the  unknown 
system  to  be  processed  separately.  This  reduces  the  sizes  of  the  data  matrices  involved 
in  the  modeling  process.  Both  block  and  sequential  processing  cases  are  formulated  but 
only  the  block  processing  case  is  tested. 

C.  OVERVIEW  OF  THESIS 

Chapter  2  is  about  ARMA  modeling.  It  also  presents  a  detailed  derivation  of  the 
method  of  ordinary  least  squares  because  it  forms  the  basis  on  which  other  modeling 
techniques  depend. 

Chapter  3  presents  a  modified  least-squares  approach  called  the  method  of  instru¬ 
mental  variables.  It  is  attractive  due  to  its  simplicity  and  good  noise  performance. 
Chapter  3  presents  results  of  using  this  method  on  several  second  and  third-order  test 
systems. 


Chapter  4  presents  a  new  multichannel  approach  to  ARMA  modeling.  This  ap¬ 
proach  is  presented  in  block  and  sequential  processing  forms.  This  chapter  also  presents 
several  adaptations  of  the  block  form  which  improve  its  speed  of  convergence. 

Chapter  5  contains  a  summary  of  the  thesis  and  lists  topics  for  further  research. 

The  appendix  contains  the  programs  used  to  test  the  sequential  IV  algorithm  and 
the  block  multichannel  iterative  algorithm.  Subroutines  common  to  both  programs  are 
grouped  together  and  listed  at  the  end  of  the  appendix. 
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II.  ARMA  MODELING 


A.  ARMA  PROCESSES 

Modeling  as  an  autoregressive-moving-average  (ARMA)  process  has  the  potential 
for  achieving  a  close  fit  to  the  system  using  a  reduced  order  over  that  which  a  moving 
average  or  an  autoregressive  model  alone  could  achieve.  ARMA  modeling  is  concerned 
with  finding  a  set  of  AR  parameters  and  MA  parameters  which  combined  describe  an 
ARMA  process  that  approximates  the  characteristics  of  a  target  system. 

The  general  form  of  the  ARMA  model  is  shown  in  Figure  2  on  page  6.  The  output 
at  time  n,  y(n),  is  a  linear  combination  of  past  outputs  and  past  and  present  inputs. 
The  a,  and  b.  are  constants  referred  to  as  tap  weights.  The  a,  parameters  form  the  MA 
part  of  the  ARMA  model.  The  b ,  parameters  form  the  AR  part.  In  equation  form  the 
output  of  the  ARMA  system  is  represented  by  the  following  difference  equation: 

,v  m 

>•(«)  -  -  Yjb<y(n  - o + _  (2- 1 ) 

1=1  1=0 

where  X  is  the  order  of  the  AR  part  of  the  ARMA  model  and  M  is  the  order  of  the  MA 
part  of  the  ARMA  model.  This  means  the  ARMA  output  at  the  current  time  depends 
on  the  last  X  values  of  the  ARMA  output.  The  X  b,  weighting  parameters  determine 
exactly  how  the  new  output  depends  on  the  past  outputs.  The  M  a,  weighting  parame¬ 
ters  determine  how  the  new  output  depends  on  the  current  and  M  —  1  past  inputs. 

Equation  (2.1)  in  vector  form  becomes: 

y  =  xT0  (2.2) 

where  x  is  a  (X  +  M+  1)  x  1  vector  of  input  and  output  data  values  given  by: 

x  =  [  -y(n  -  1)  -y(n-2)  ...  —y(n  -  X)  x(n)  *(«-!)  ...  x{n  -  Xf)]T  ( 2.3) 

and  6  is  a  {X  +  M  +  1)  x  1  vector  of  the  AR  and  MA  tap  weights  given  by: 

6  =  [6,  b2  ...  btf  Oq  ay  ...  cw]  (2-4) 


5 


Figure  2.  Structure  of  the  ARM  A  model 


For  N  +  L  —  1  data  points  available  for  y  and  M  +  L  data  points  available  for  u,  we 
can  write  a  block  equation  which  gives  the  value  of  the  output  at  progressive  sampling 
times: 


*1 

*2 


y(n-L+  1) 

xT{n  +  1) 

An  -  L  +  2) 

xT(n  +  2) 

xr(n  +  L) 

(2.5) 


The  /*"  row  in  equation  (2.5)  is  the  value  ofy  at  time  n  +  i  based  on  output  data  available 
through  time  n  +  /  —  1  and  input  data  available  through  n  +  /.  The  <*h  row  is  identically 
equation  (2.2)  at  time  n  +  i.  In  vector  form  equation  (2.5)  becomes: 

y*X0  (2.6) 

where  0  is  defined  in  equation  (2.4);  y,  the  vector  of  output  values,  is  given  by: 

y  =  Ot«-L+l)  A"-L  +  2)  ...  An)f  (2.7) 


and  X  is  a  partitioned  matrix  with  rows  comprised  of  data  vectors  exactly  like  equation 
(2.3)  only  shifted  in  time.  At  successive  sampling  times,  when  new  data  is  obtained,  data 
used  to  calculate  the  previous  output  shifts  one  column  to  the  right.  The  new  data  fills 
in  the  left  mosty  and  u  columns.  The  matrix  X  is  given  by: 


-yin  ~  L) 

—y{n  —  L+  1) 


I 


-yin-  l) 


—An~  +1)  «(«—£+  I) 
~An  —  fl+  2)  u(n  —  L  +  2) 


-An  ~  N)  u(n ) 


k(h-m+1) 
u(n-  n  +2) 


(2-8) 


u{n  —  M) 


where  rj  is  defined  as  A'  +  L  and  n  is  defined  as  M  +  L. 
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If  the  a,  and  b,  are  estimates  of  the  true  values  of  the  AR  and  MA  parameters,  then 
the  filter  output  will  be  an  estimate  of  the  true  output.  We  use  a  hat  over  a  variable  (for 
example,  y)  to  indicate  an  estimated  value.  Rewriting  equation  (2.6)  using  the  estimated 
ARMA  parameters  results  in: 

y  =  Xe  (2.9) 

A 

where  6  is  defined  as: 


e  =  [6,  b2  ...  bN  3o  «1  ...  (2.10) 

and  y  is  now  the  vector  of  estimated  output  values  and  is  given  by: 

}'  “  CP(«  —  L  +  l)  y(n  —  L  + 2)  ...  y(«)]r  (2.11) 

L'p  until  now,  we  have  discussed  estimating  the  output  of  a  system  given  its  input, 
past  output,  and  an  estimate  of  the  parameters  which  describe  it.  If,  however,  we  know 
the  output  and  input  of  the  system,  based  on  these  equations,  we  can  use  them  to  gen- 

a 

erate  a  set  of  a,  and  6,  which  produces  an  ARMA  output  which  is  the  best  possible  es- 

A  A 

timate  of  the  system  output.  Then  the  a,  and  b,  wrill  be  optimal  parameters  for  describing 
the  operation  of  the  unknown  system  as  an  ARMA  process. 

B.  METHOD  OF  ORDINARY  LEAST-SQUARES 

In  this  thesis  we  use  the  method  of  ordinary  least-  squares  as  the  means  of  finding 
the  optimum  set  of  ARMA  parameters.  It  is  a  well  known  modeling  technique.  It  offers 
the  advantage  of  being  widely  used  in  the  scientific  community  for  a  variety  of  modeling 
problems.  It  has  been  applied  successfully  to  a  large  number  of  modeling  problems  with 
good  results  and  has  been  successfully  applied  to  classes  of  problems  for  which  other 
methods  have  failed.  (Ref.  3:  p.  4) 

To  apply  the  method  of  ordinary  least-squares  to  system  identification  we  form  the 
error  between  the  actual  system  output  and  the  estimated  output  generated  by  the 
ARMA  model.  This  error  is  given  by: 

«-y-y-y-»  (2.12) 

where  y  is  the  vector  of  the  actual  system  outputs  given  by  equation  (2.7)  and  y  is  the 
vector  of  ARMA  outputs  given  by  equation  (2.11).  [Ref.  1:  p.  176] 
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We  then  let  the  sum  of  the  squares  of  the  errors  at  the  instances  of  time  the  meas¬ 
urements  of  the  data  were  taken  become  a  measure  of  how  well  the  estimates  approxi¬ 
mate  the  true  system  outputs.  This  measure  of  performance,  or  cost  function,  is  denoted 
J  .  It  is  written  in  equation  form  as: 


n+L 


lan+l 


Replacing  the  error  in  equation  (2.13)  with  its  equivalent  expression  from  equation 
(2.12)  results  in: 

T"  At*  tv*‘  A  t* 

J=  yry  +  0rXXr0  -  20rXy  (2.14) 

Equation  (2.14)  shows  that  the  performance  measure  is  a  function  of  the  estimated  pa¬ 
rameters.  The  criterion  is  to  minimize  the  measure  of  performance  by  taking  its  deriva¬ 
tive  with  respect  to  the  parameter  estimates  and  setting  it  equal  to  zero.  Then  equation 
(2.14)  becomes: 

-  0  =  0  +  2XX70  -  2Xry  (2.15) 

60 

A 

Solving  for  0,  the  parameters,  gives  us  the  result: 

0  =  (XTX)“,X7y  (2.16) 

Equation  (2.16)  is  the  ordinary  least-squares  solution  for  the  optimum  ARM  A  pa¬ 
rameters.  It  provides  the  best  possible  description,  in  a  least-squares  sense,  of  the  data 
source.  The  resulting  parameters  provide  the  closest  fit  to  the  actual  input  and  output 
data  of  the  system  in  the  sense  of  least-squares  errors. 

Equation  (2.16)  uses  a  block  processing  approach.  The  product  of  XfX  must  be 

A 

formed  and  then  inverted  in  order  to  calculate  0  .  In  addition  to  being  computationally 
intensive,  the  estimate  cannot  be  updated  when  new  data  becomes  available  without  re¬ 
calculating  (X^)-1  .  A  sequential  update  which  does  not  require  (X^)-1  to  be  recalcu¬ 
lated  is  presented  in  the  next  chapter  in  the  context  of  the  instrumental  variable  method 
of  least-squares. 
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III.  INSTRUMENTAL  VARIABLE  METHOD  OF  SYSTEM 

IDENTIFICATION 


A.  INTRODUCTION 

The  instrumental  variable  (IV)  method  of  system  identification  is  a  variation  of  the 
method  of  ordinary  least-squares.  Its  attraction  over  ordinary  least-squares  is  that  there 
is  no  bias  in  estimating  the  parameters  when  dealing  with  noise  [Ref.  4:  p.  406].  Also, 
this  method  is  known  to  yield  consistent  estimates  and  remains  as  easy  to  use  as  the 
method  of  ordinary  least-squares  [Ref.  3:  p.  119). 

When  an  additive  noise  term  is  present  in  the  observable  output,  >-(w),  the  output  is 
given  by: 

y{n)  =  w{n)  +  v(«)  (3.1) 

Here  w(«)  represents  the  actual  output  of  the  system  and  v(«)  represents  the  noise.  When 
this  noise  has  a  non-zero  mean,  using  the  noise  corrupted  output  to  model  the  unknown 
system  by  the  ordinary  least-squares  approach  leads  to  inaccurate  estimates  of  its  pa¬ 
rameters.  The  parameters  are  referred  to  as  biased  estimates.  [Ref.  3:  p.  119,  Ref.  1:  pp. 
192-193,  and  Ref.  5:  p.  704). 

The  IV  method  shown  in  Figure  3  on  page  1 1  generates  an  estimate  of  the  un¬ 
known  system's  output  by  processing  the  input  data  through  an  auxiliary  model  which 
closely  approximates  the  unknown  system.  In  our  implementation  of  the  IV  method, 
the  auxiliary  model  is  an  ARMA  model.  Its  output  is  free  of  the  noise  affecting  the 
unknown  system.  The  IV  method  uses  the  auxiliary  model  output  (estimate),  w,  to  cal¬ 
culate  the  parameters  of  the  unknown  system,  Therefore  the  IV  parameter  estimates  are 
not  biased  like  those  generated  by  the  method  of  ordinary  least-squares. 

The  IV  method  assumes  the  existence  of  a  matrix  Z  composed  of  the  auxiliary 
model's  input  and  output  data  which  has  the  following  two  properties  [Ref.  4:  p.  406): 

Jim  “aT  Z7*  =  0  (3.2) 

A-*oo  A 

lim-^ZrX  =  Q  (3.3) 

\-*00  4  ’ 

where  t  is  the  error  in  fitting  the  parameter  estimates  to  the  data  and  is  given  by: 
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v(n) 


Figure  3.  Modeling  by  the  instrumental  variable  method 

«  *  y  -  Xd/y  (3.4) 


and  Q  is  a  nonsingular  square  matrix. 

The  First  property  means  Z  is  orthogonal  to  the  error.  This  leads  to  the  cancellation 
of  the  bias  term  inherent  in  ordinary  least-squares  techniques  (Ref.  4:  p.  406].  The  sec¬ 
ond  property  ensures  the  inverse  of  ZK  exists.  Z  is  assumed  to  have  the  same  structure 
and  size  as  the  data  matrix  X  in  equation  (2.8).  Its  contents  dilfer  in  that  the  noise 
corrupted  system  output  >•(«)  in  X  is  replaced  by  the  output  of  the  auxiliary  model  w(n) 
in  Z.  The  new  data  matrix  Z  is  given  by: 


—  w(n  —  L) 

-  w{n  -  L  +  I) 


-  w(n  -  r\  +1) 

—  w(n  -  7  +2) 


u(n  —  L  +  1) 
u(n  —  L +  2) 


u(n  —  n  +1) 
u{n  —  n  +2) 


(3.5) 


—  \v{n  —  1)  ...  —  w(n  —  A")  u(n) 


u{n  —  M) 
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where  tj  is  defined  as  S  +  L  and  n  is  defined  as  Af  +  L.  Comparing  X  in  equation  (2.8) 
and  Z  in  equation  (3.5),  we  note  the  substitution  of  w{n)  for  >•(«).  Thus,  we  are  now 
using  estimates  of  the  true  output  \v(n)  instead  of  noise  corrupted  samples  >•(/?). 

To  incorporate  Z  into  the  parameter  estimation  process  we  begin  with  equation 
(2.12),  which  we  rewrite  as: 

y  =  XO  +  e  (3.6) 

A 

This  equation  says  that  the  estimates  of  the  output,  given  by  X0,  differ  from  the  actual 
outputs,  y,  by  some  fitting  error  «.  Multiplying  equation  (3.6)  by  Zr  yields: 

ZTy  =  ZrX0  +  Z  Tt  (3.7) 

Equation  (3.3)  ensures  that  Z*X  can  be  inverted.  Solving  for  0  results  in: 

6  =  (ZrX)-]Zry  -  (ZrX)~}Zre  (3.8) 

The  (Z^-'Z7}'  term  in  equation  (3.8)  is  the  IV  estimate  of  the  parameters.  It  is  written 
as: 

0/v  =  (ZrXr'Zry  (3.9) 

The  (Z^-'Z^  term  in  equation  (3.8)  represents  a  potential  bias  in  the  estimate.  The 
first  property  of  the  Z  matrix,  given  in  equation  (3.2),  ensures  this  bias  goes  to  zero, 
asymptotically.  Applying  this  property,  equation  (3.8)  can  be  rewritten  as: 

6  =  (ZTX)-'Zry  =  6ly  (3.10) 

Equation  (3.10)  gives  an  unbiased  estimate  of  the  ARMA  parameters.  [Ref.  1:  pp. 
192-193]: 

Other  least-squares  methods  avoid  the  bias  inherent  in  ordinary  least-squares  but 
they  are  more  complicated  than  the  IV  method  to  implement  [Ref.  3:  p.  119].  Although 
this  thesis  does  not  attempt  an  analysis  of  the  IV  method  in  the  presence  of  noise,  any 
practical  system  identification  technique  must  deal  with  noise.  Hence,  the  attraction  of 
and  the  desire  to  use  the  IV  method. 

Equation  (3.10)  represents  the  block  processing  case.  It  assumes  N  4-  L  —  1  output 
samples  and  M  +  L  input  samples  are  available.  These  samples  are  used  to  calculate  an 


estimate  of  the  parameters.  Samples  beyond  this  range  are  not  included  in  the  esti¬ 
mation  process.  Block  processing  involves  multiplication  of  two  L  x  (A'+  M  +  1)  ma¬ 
trices  to  form  a  third  matrix.  Then  this  third  matrix  must  be  inverted.  This  is  a 
computationally  intensive  process.  In  what  follows,  we  present  a  sequential  algorithm 
to  compute  6IV  which  avoids  matrix  inversions. 

B.  SEQUENTIAL  LEAST-SQUARES  ESTIMATION  USING  INSTRUMENTAL 
VARIABLES 

A  sequential  process  for  estimating  the  parameters  of  an  unknown  system  requires 
fewer  computations  than  a  block  process.  In  a  manner  similar  to  that  presented  in  Hsia 
[Ref.  3:  pp.  22-25]  for  the  general  least-squares  case,  the  block  IV  estimation  process 
described  above  can  be  converted  into  a  sequential  IV  estimation  process.  Using  the 
sequential  process  also  allows  the  coefficients  to  be  updated  based  on  the  new  data  that 
becomes  available. 

The  derivation  of  the  sequential  estimation  procedure  consists  of  two  parts.  The 
first  part  is  the  derivation  of  an  equation  to  update  the  data  matrix,  Q {m  +  1) ,  based 
on  the  previous  data  matrix,  Q(m) ,  and  the  new  data:  w{m) ,  j’fm) ,  and  u{m  +  1) 
where  m  represents  the  iteration.  The  second  part  involves  developing  an  equation  for 
updating  the  estimate  of  the  parameters,  0/v(m  +  1),  based  on  the  previous  estimate, 
0n{m) ,  the  previous  data  matrix,  Q(m) ,  and  the  new  data:  w(m) ,  y{m) ,  and 
u{m  +  1)  . 

Define  the  data  matrix  Q (m)  to  be: 

QW  =  [Zmrxm]-’  (3.11) 

where  Zm  is  given  by  equation  (3.5)  and  X„  is  given  by  equation  (2.8).  The  property  of 
equation  (3.3)  assures  that  Q  exists.  Note  that  Q (m)  includes  output  data  available 
through  m  and  input  data  available  through  m  +  1  .  Since  both  Z„  and  X„  are 
m  x  (A’  +  M)  matrices,  Q(m)  will  be  a  (Ar  +  M)  x  (N  +  Af)  matrix.  As  the  number  of  rows 
of  Z  and  X  increase  to  accommodate  the  increasing  numbers  of  data  points,  the  size  of 
Q  will  remain  the  same.  At  the  next  sample  time,  i.e.,  at  m  +  1,  the  data  matrix  becomes: 

Q(m+l)-[Zj+,X„+,r'  (3.12) 

where  the  data  matrices  at  m  +  l  are  given  by: 
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I 


I 

p 


^■m+ 1  — 


zT(m  +  1) 


(3-13) 


( 

I 


x„„- 


xT{m  +  1) 


(3.14) 


and  z T(m+  1)  and  xT{m  +  1)  are  vectors  which  contain  the  most  recent  data  values. 
Substituting  equations  (3.13)  and  (3.14)  into  equation  (3.12)  and  expanding,  results  in: 

(3-15) 

Expanding  further  yields: 

Q(m  +  1)  -  [ZX  +  z(m  +  l)xr(m  +  I)]-1  (3.16) 

In  equation  (3.16)  we  see  that  two  terms  make  up  the  new  data  matrix.  The  Z£Xm  term 
is  all  the  data  that  was  available  through  time  m.  The  z (nt  +  l)xr(m  +  1)  term  contains 
the  new  data.  To  perform  the  inversion,  let  A  =  Zptn,  B  =  z(m+1),  C=1  and 
D  =  xT(m  +  1)  .  Then  by  the  matrix  inversion  lemma: 

Q(m  +  1)  -  A-'  -  A-1B(C-1  +  DA",Br’DA'1  (3.17) 

Substituting  the  appropriate  expressions  for  A,  B,  C,  and  D  back  into  equation  (3.17) 
yields  the  equation: 

Q(m  +  1)  =  (Z  X)'1  -  (ZX)''^  +  1) 

.  [l+xV+lXZXrVm  +  l)]-' 

.  xT(m  +  l)(ZXr‘ 

Substituting  Q(m)  for  (Z£XJ-'  reduces  equation  (3.18)  to: 

Q(m  +  1)  =  Q(m)  -  Q(m)z(m  +  1)[1  +  xT(m  +  l)Q(m)z(m  +  1)]“’ 

•  xr(m+l)Q(m) 


(3.18) 


(3.19) 


Q(«+  l)  = 


X  z("J+l)] 


X 


m 


xT{m  +  1) 
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This  completes  the  first  step  of  the  derivation.  Equation  (3.19)  expresses  Q  at  time 
m- 1- 1  in  terms  of  the  old  Q  and  the  new  data.  The  term  in  the  brackets  is  a  scalar. 
Computational  intensity  has  been  reduced  because  a  large  matrix  does  not  have  to  be 
generated  and  its  inverse  does  not  have  to  be  calculated. 

A 

Continuing  with  the  derivation,  the  estimate  6IV  for  data  available  through  m  can 
be  written  as: 


Mm)  =  (ZlXmr’zIym  (3.20) 

The  estimate  6IV  for  data  available  through  m  +  1  can  be  written  as: 

e,M  +  1)  -  (Z£+!Xm4,r’z£+Iym+I  (3.21) 

Substituting  equation  (3.12)  into  equation  (3.21)  results  in  an  expression  for  the  estimate 
of  the  parameters  in  terms  of  the  new  data  matrix  and  all  the  available  data  given  by: 

+  1)  =  Q(m  4-  l)Z£+1yOT+1  (3.22) 

M,-i+l)  =  Q(m+l)[z£  z{m  +  1)]  -  (3.23) 

y{m  +  1) 

+  1)  =  Q(m  +  l)[Z£ym  +  z  (m  +  1M"»  +  1)3  (3-24) 

Substituting  for  Q(m  +  1)  from  equation  (3.19)  and  expanding  results  in: 

eiy(m  +  l)  =  Q(m)z£ym 


-  Q(m)z(m  +  1)[1  +  xT(m  +  l)Q(m)z(m  +  l)]'!jrr(m  +  l)Q(m)Z£ym 
+  Q(m)z(m  +  l)y(m  +  0  (3-25) 

—  Q(m)z(m  +  1)[1  +  xT(m  +  1  )Q(m)z(m  +  1)]_1 
•  xT(m  +  1)Q  (w)z(m  +  1Mm  +  1) 

Although  somewhat  lengthy,  this  equation  has  the  desired  form.  To  simplify  it,  its  last 
two  terms  can  be  arranged  into  the  form: 

Q(m)z(m  +  1){1  -  [1  +  xr(m  +  I)Q(m)z(m  +  l)]-1xr(m  +  l)Q(m)z(m  4-  1)}  .- 
•  y{m  +  1) 
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The  terms  within  the  braces  can  be  thought  of  as  the  result  of  a  previous  application  of 
the  matrix  inversion  lemma  with  A-‘  =  l,  B=I,  C_l  =  1  and 
D  =  \r{m  +  l)Q(m)z(wt  +  1)  .  Reversing  the  lemma  results  in: 

Q(m)z(m  4-  1)[1  +  \T{m  +  l)Q(m)z{m  -f  +  1)  (3-27) 

Replacing  the  last  two  terms  in  equation  (3.25)  with  this  result  gives  us: 

6iy(m  +  1)  =  Q(m) Z^ym  -  Q(m)z(m  +  1)[1  +  xr(m  +  l)Q(m)z(/n  +  l)]"1 

.  x\m  +  l)Q(m)Z£ym  +  Q(m)z(m  +  1)  (3.28) 

•  [1  +  xT{m  +  l)Q(m)z(m  +  1)]“  +  1) 

Factoring  Q(m)z(m  +  l)  and  [1  +  xT{m  +  l)Q(m)z(m  4-  1)]_1  from  the  last  two  terms  re¬ 
duces  equation  (3.28)  further  to: 

6n(m  +  1)  -  Q(m) Z£ym  +  Q(m)z(m  +  l)[l  +  xT(m  +  l)Q(m)z(m  +  1)]“'  29 

•  [y\m  +  1)  ~  +  i)Q(«)Z*ym] 

Substituting  equation  (3.11)  into  equation  (3.20)  and  then  equation  (3.20)  into  equation 
(3.29)  yields  the  final  form  for  the  update  of  the  estimate  of  the  parameters: 


dlv{m  4-  1)  =  On{rn)  4-  Q(m)z(m  +  1) 

•  [14-  xr(m  4-  l)Q(m)z(m  +  I )]~ '  CKm  4-  1)  —  xT(m  4-  1  )0/(X*t)] 


(3.30) 


This  is  the  desired  result  for  updating  the  estimate  of  the  parameters.  Note  that  like 
equation  (3.19),  the  matrix  inversion  of  equation  (3.21)  has  been  reduced  to  inversion 

A 

of  a  scalar.  Equation  (3.30)  describes  the  update  of  6/v(ni  4-  1)  in  terms  of  the  previous 

A 

estimate  of  the  parameters,  the  previous  data  matrix,  Q(m),  and  the  new  data: 

w(m) ,  y(m) ,  and  u{m  4-  1) . 

C.  TESTING  THE  SEQUENTIAL  INSTRUMENTAL  VARIABLE  ALGORITHM 
Equations  (3.19)  and  (3.30)  above  comprise  the  sequential  IV  algorithm.  Several 
tests  of  this  algorithm  were  made  using  second  and  third-order  filters  as  unknown  sys¬ 
tems.  Tests  were  run  via  computer  simulation  using  the  filters  to  generate  the  output 
data.  A  Gaussian  random  process  with  zero  mean  and  unit  variance  was  used  as  the 
input.  The  input  was  produced  by  1MSL  subroutine  GGNML.  Graphs  were  created 
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using  DISSPLA.  Table  1  on  page  17  shows  pole  and  zero  locations  as  well  as  numera¬ 
tor  and  denominator  parameters  for  the  test  filters. 


Table  1.  TEST  SYSTEMS  FOR  THE  IV  MODELING  METHOD 


TEST 

LOCATIONS 

LOCATIONS 

AR 

MA 

FILTER 

OF  POLES 

OF  ZEROS 

PARAMETERS 

PARAMETERS 

T2 

0.445  + jO.228 

0.4  + j  1.273 

1.0 

0.5 

0.445  -  jO.228 

0.4-  j  1.273 

-0.89 

-0.4 

0.25 

0.S9 

T2N 

0.445  + jO.228 

0.4  + jO. 8 

1.0 

1.0 

0.445  -jO.228 

0.4  -j0.8 

-0.89 

-0.80 

0.25 

0.80 

T3 

0.6605 

-1.0 

1.0 

0.6647  +  j0.502 

-1.0 

-1.99 

0.6647-j0.502 

-1.0 

1.57 

-0.45S 

0.0154  i 

Results  of  the  tests  are  shown  in  graphical  form  in  Figure  4  on  page  18,  Figure  5 
on  page  19,  and  Figure  6  on  page  20.  Dashed  lines  indicate  the  true  values  of  the  pa¬ 
rameters.  Solid  lines  are  the  IV  method's  estimates. 

For  both  second-order  test  cases  shown,  the  algorithm  converged  quickly  and 
produced  accurate  results.  For  the  third-order  test  case,  convergence  took  longer  but 
the  values  were  accurate.  A  third-order  system  is  more  complex  than  a  second-order 
system,  so  conceivably  it  would  require  more  iterations  to  converge.  The  number  of  it¬ 
erations  required  is  of  the  same  order  as  the  method  of  ordinary  least-squares. 

Table  2  on  page  21  contains  the  IV  algorithm  s  best  estimates  of  the  parameters  and 
the  number  of  iterations  required  to  converge  to  those  estimates.  It  also  shows  the  ab¬ 
solute  and  percent  errors  from  the  true  parameters. 


PARAMETER  VALUES 


ITERATIONS 


(A) 


ITERATIONS 


(B) 

Figure  6.  Tliird-order  test  case  T3.  (A)  MA  parameters.  (B)  AR  parameters. 
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Tahle  2.  COEFFICIENT  ESTIMATES  BY  THE  IV  MODELING  METHOD. 


IV.  SYSTEM  IDENTIFICATION  USING  AN  ITERATIVE 
MULTICHANNEL  APPROACH 

A.  INTRODUCTION 

This  chapter  presents  an  alternate  system  identification  method,  the  iterative  multi¬ 
channel  approach.  This  approach  differs  from  the  IV  method  and  the  method  of  ordi¬ 
nary'  least-squares  presented  in  the  preceding  chapters  in  that  it  processes  the  input  and 
output  data  from  the  unknown  system  in  separate  channels.  In  its  block  processing 
form  one  advantage  over  the  IV  and  ordinary  least-squares  methods  is  a  reduction  in  the 
sizes  of  the  data  matrices.  As  a  result,  the  computational  complexity  of  the  multichannel 
algorithm  is  on  the  order  of  M2  +  A"2  .  where  M  is  the  order  of  the  MA  part  and  A'  is  the 
order  of  the  AR  part.  In  contrast,  the  block  IV  and  ordinary  least-squares  methods  re¬ 
quire  computations  on  the  order  of  (A/  +  N)1 . 

B.  PREVIOUS  MULTICHANNEL  METHODS 

Whittle  [Ref.  6:  pp.  129-130]  was  the  first  to  develop  a  multichannel  solution  for  the 
ARMA  modeling  problem.  He  sought  to  extend  the  recursive  Durbin  solution  for  esti¬ 
mating  the  parameters  of  a  single  variable  autoregressive  process  to  a  multivariable 
autoregressive  process.  He  discovered  that  to  do  this  he  would  have  to  fit  the  data  to 
two  autoregressive  processes  simultaneously.  One  of  the  autoregressions  would  use 
present  data  samples  to  predict  the  value  of  the  data  one  time  step  in  the  future.  This 
is  called  forward  prediction.  The  second  autoregression  would  use  present  data  samples 
to  predict  the  value  of  the  data  at  the  previous  time  instant  and  is  referred  to  as  back¬ 
ward  prediction.  Sometime  during  this  research,  Whittle  determined  that  if  the  input 
was  derived  from  a  MA  scheme,  making  the  process  ARMA,  then  the  solution  w'ould 
remain  the  same  provided  the  correlations  of  the  input  used  in  the  parameter  estimation 
process  had  shifts  greater  than  the  MA  scheme.  Whittle's  use  of  the  two  separate  and 
simultaneous  autoregressions  to  model  an  ARMA  process  can  be  thought  of  as  a 
multichannel  modeling  approach. 

Further  work  in  the  area  of  multichannel  ARMA  modeling  was  conducted  by  Perry 
and  Parker  [Ref.  7:  pp.  509-510],  They  started  out  with  the  ARMA  problem  formulation 
discussed  in  Chapter  2.  Using  the  method  of  ordinary  least-squares  to  minimize  the 
mean  square  error,  they  found  the  solution  for  the  estimate  of  the  parameters  to  be  the 
Wiener  solution  given  by: 
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(4.1) 


Ryu'  lm  K 

Ry'y  RyvJLa  J  Lr>-' 

In  equation  (4.1)  R„  is  a  matrix  of  autocorrelations  of  the  past  outputs,  R„.„.  is  a  matrix 
of  autocorrelations  of  the  inputs,  Ryu,  and  R^  are  crosscorrelations  of  the  input  and 
output  data,  b  is  a  vector  of  AR  parameters,  and  a'  is  a  vector  of  MA  parameters.  In 
addition,  r„  is  a  vector  of  autocorrelations  of  past  output  data  with  the  current  output 
and  r>u.  is  a  vector  of  crosscorrelations  of  input  data  with  the  current  output.  By  as¬ 
suming  the  first  MA  parameter,  a\  was  known,  they  were  able  to  treat  it  as  a  gain  and 
factor  it  out  of  all  the  other  MA  parameters.  This  allowed  them  [Ref.  7:  pp.  509-510] 
to  extract  the  (;V+  1)"  row  and  column  of  equation  (4.1)  and  rewrite  the  solution  in  the 
form: 

rRyy  Mm  m  ivi  (4.2) 

j^RUy  RUUJL«J  L*VuJ  LruuJ 

In  equation  (4.2),  a'  has  been  rewritten  as  a  to  indicate  that  the  a'0  term  has  been  fac¬ 
tored  out  of  all  of  the  MA  parameters.  Reasoning  that  equation  (4.2),  the  ARM  A  sol¬ 
ution,  was  a  generalization  of  the  AR  solution,  they  figured  that  it  must  have  a  recursive 
solution  consisting  of  some  combination  of  the  Levinson-Durbin  algorithm,  a  recursive 
solution  for  the  AR  problem.  They  then  determined  equation  (4.2)  was  in  a  form  similar 
to  Whittle's  formulation  of  the  problem.  So  they  reasoned  that  they  could  use  a  form 
of  Whittle's  solution  to  solve  the  ARM  A  modeling  problem.  Like  Whittle,  their  solution 
consists  of  a  forward  and  a  backward  autoregression.  It  uses  two  coupled  lattice  filters 
to  process  the  input  and  output  data.  OIT  diagonal  elements  of  the  lattice  coefficient 
matrices  specify  the  coupling  points  of  the  two  lattices. 

C.  ITERATIVE  APPROACH  TO  MULTICHANNEL  ARMA  MODELING 

This  thesis  proposes  another  solution  to  the  ARMA  modeling  problem  using  the 
multichannel  approach.  It  is  an  iterative  approach  with  no  direct  coupling  of  the  two 
channels.  However,  note  that  there  is  an  implicit  coupling  in  the  sense  that  the  ARMA 
system  output  samples  yfn)  are  a  function  of  the  present  and  past  input  samples  u(n), 
and  the  past  output  samples.  This  is  shown  in  equation  (2.1).  The  approach  proposed 
uses  two  separate  fmite-impulse-response  (FIR)  filters  to  estimate  the  unknown  system. 
One  filter  estimates  the  AR  part  of  the  unknown  system.  The  second  estimates  the  MA 
part  of  the  unknown  system.  Figure  7  on  page  25  shows  the  structure  of  this  approach. 
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The  derivation  of  the  equations  for  this  approach  follows  the  method  of  ordinary  least- 
squares. 

A(z) 

From  Figure  7  on  page  25,  the  value  of  the  signal  F(z)  is  seen  to  be  — — -  l'(z)  . 

»(-) 

When  this  signal  passes  through  C(z) ,  if  C(z)  is  close  to  B(z),  the  resulting  signal  A’(z)  is 
approximately  A (z)L'(z).  Also  from  Figure  7  on  page  25,  the  value  of  Z(z)  is  seen  to  be 
A(z)U(z)  provided  D(z)  is  dose  to  A (z).  The  difference  of  these  two  signals  forms  the 
error  which  we  seek  to  minimize  by  the  method  of  ordinary  least-squares.  In  minimizing 
the  error  we  seek  to  drive  D(z)  and  C(z)  as  close  as  possible  to  A(z)  and  B(z)  ,  respec¬ 
tively. 

Referring  to  Figure  7  on  page  25,  signals  z{n)  and  x(x)  are  defined  as  the  outputs 


from  two  FIR  filters  and  are  given  by  the  equations: 

z(n)  =  d0u(n)  +  d\u(n  —  I)  +  d2u(n  —  2)  +  —  +  dMu(n  —  Af)  =  ur(«)d  (4.3) 

x(n) -y{ n)  +  -  1)  +  c-#{n  -  2)  +  -  +  -  N)  =  yT{n)c  (4.4) 

where  the  vectors  d  and  c  represent  the  systems  D(z)  and  C(z),  respectively.  The  vectors 
d,  u(w),  c,  and  y(n)  are  given  by  the  following  equations: 

d  -M,  4  d2  ...  dMf  (4.5) 

u(w)  =  [«(«)  «(«-])  «(« —  2)  ...  (4.6) 

C  —  Cl  Cl  c2  ...  c*]7'  (4.7) 

>’(«)  =  W»)  ><«-•)  -  2)  ...y(n  -  A^]r  (4.8) 


The  d  parameters  are  estimates  of  the  MA  portion  of  the  ARMA  process.  The  c  pa¬ 
rameters  approximate  its  AR  portion.  The  vector  u(n)  is  the  input  data  vector  of  length 
M,  the  order  of  the  MA  part,  and  y(n)  is  the  output  data  vector  equal  of  length  N,  the 
order  of  the  AR  part. 

Forming  the  error  between  x  and  z  results  in: 

e(n)  —  z(n)  —  x (n)  =  ur(n)d  —  yr(«)c  (4.9) 

The  least-squares  performance  criterion  is: 
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UNKNOWN  SYSTEM 


e(n)  -  z(n)  -  x(n) 


Figure  7.  Multichannel  modeling  process 


L  L 


xin)f 


(4.10) 


>1*1  nm  | 

Substituting  for  e{n )  from  equation  (4.9),  we  can  write  equation  (4. 10)  in  vector  form  as: 

L 

y=^(uT d-yV  (4.11) 


where  we  have  dropped  the  time  index,  n,  for  convenience.  Expanding  equation  (4.11) 
results  in: 


L 

^druurd  +  cTyyTc  -  2druyrc 

n*1 


(4.12) 
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We  notice  that  the  performance  criterion  is  a  function  of  both  d  and  c.  Minimizing  the 
performance  criterion  by  differentiating  with  respect  to  the  vector  c  and  equating  the 
results  to  zero  yields: 


L  L 

-§f  =  0  =  0  + 2^T(yy7)c -  2^Myur)d 


(4.13) 


/»=  1  ,  /!=) 

Similarly,  differentiating  the  performance  criterion  with  respect  to  the  vector  d  and 
equating  the  result  to  zero  yields: 


8J_ 

c  d 


L  L 

=  0  =  0  +  2^T(uur)d  -  2^(uyr)c 

7i«l  n«l 


(4.14) 


Solving  equation  (4.13)  for  c  and  equation  (4.14)  for  d  results  in  two  equations  for 
estimating  the  AR  and  MA  parameters  of  the  unknown  system  given  by: 

c  =  R“'R^d  (4.15) 


and 


where 


d  =  R  J  Ru>c 


n=l 

ruu(0)  ruu(  1)  .  .  ruu(i\) 
ruuW  ruu(°)  ■  ■ 

rn{N)  .  .  .  rm{ 0) 

(4.16) 


(4.17) 


is  an  estimate  of  the  input  autocorrelation  matrix.  The  elements  of  R„  are  computed 
using  the  unbiased  method  as  follows: 

L—l 

=  7^  £*/(/>(/' -o  for  /  =  0,  1,  2,  ...  ,M  (4.18) 

o 
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Matrices  Ryy,  R„y,  and  Ryv  appearing  in  equations  (4.15)  and  (4.16)  have  structures  iden¬ 
tical  to  equation  (4.17),  where  Ryy  is  the  estimate  of  the  output  autocorrelation  matrix, 
and  Riy  and  Ryk  are  estimates  of  the  crosscorrelation  matrices.  Note  that  Ryu  =  R£  .  The 
elements  of  these  matrices;  ryy,  r„y  and  ry„,  are  computed  as  follows: 

L~l 

ryyW  =  J^Y/jW-0  '-0,1,  2,  ...  ,  N  (4.19) 

y-o 

L-t 

for  /  =  0,  1,  2,  ...  ,  M  (4.20) 

J~ o 

and 

L-l 

VW-T for  /-0,  1,  2,  ...  (4.21) 

y=o 

Up  to  this  point,  following  the  standard  least-squares  procedure  has  resulted  in  two 
dependent  or  coupled  equations  to  solve  for  the  parameters  of  an  unknown  system 
modeled  as  an  ARMA  process.  How  best  to  solve  these  equations?  By  iteration.  The 
steps  of  the  iterative  process  are  to 

•  Calculate  the  correlation  matrices  and  vectors  from  the  available  data. 

•  Initialize  c  .  Exploit  the  fact  that  the  first  parameter  in  c  is,  c„  =  1. 

•  Solve  for  d  from  this  initial  c  . 

•  Solve  for  c  from  d  . 

•  Repeat  beginning  at  the  third  step. 

Here  is  a  summary  of  the  equations  in  proper  order  for  implementing  the  algorithm: 

•  Compute  R„, ,  Ryy  and  R,y  from  equations  (4.17)  to  (4.21).  Note  that  R,,  =  R£. 

•  Initialization: 

c(0)  =  [l  0  ...  Of  (4.22) 

•  For  k  =  0  to  K 

d(*+1)  =  R^'V^  (4.23) 

c(*+I)  =  R"IRjmd(*+l)  (4.24) 

where  k  is  the  iteration  number. 
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This  is  a  very  simple  algorithm.  For  the  case  where  the  AR  and  MA  orders  are  equal, 
the  correlation  matrices  are  half  the  size  of  the  block  data  matrices  which  must  be  gen¬ 
erated  and  inverted  in  the  IV  algorithm. 

1.  Testing  the  Multichannel  Iterative  Algorithm 

We  tested  the  algorithm  by  computer  simulation  using  second  and  third-order 
filters  as  unknown  systems.  Table  3  shows  pole  and  zero  locations  as  well  as  MA  and 
AR  parameters  for  the  test  filters.  Data  lengths  of  500  data  points  were  used  to  calculate 
the  autocorrelation  and  crosscorrelation  matrices.  Besides  the  reported  cases,  we  tested 
the  algorithm  on  several  other  second,  third  and  mixed-order  cases.  The  performance 
of  the  algorithm  in  all  cases  that  we  tested  was  fairly  uniform. 


Table  3.  TEST  SYSTEMS  FOR  ITERATIVE  MULTICHANNEL  MODELING 
METHOD 


TEST 

LOCATION 

LOCATION 

AR 

MA 

FILTER 

OF  POLES 

OF  ZEROS 

PARAMETERS 

PARAMETERS 

T2 

0.445  + jO.228 

0.4  4-j  1.273 

1.0 

0.5 

0.445  -  j0.22S 

0.4  -j  1.273 

-0.S9 

-0.4 

0.25 

0.S9 

T2X 

0.445 +  j0.228 

0.4  + j0.8 

1.0 

1.0 

0.445  -j0.228 

0.4  -j0.8 

-0.89 

-0.80 

0.25 

o.so 

T3 

0.6605 

-1.0 

1.0 

0.0154 

0.6647 +  j0. 5020 

-1.0 

-1.99 

0.0462 

0.6647  -jU.5020 

-1.0 

1.572 

0.0462 

-0.45S3 

0.0154 

Results  of  the  tests  are  shown  in  graphical  form  in  Figure  S  on  page  29, 
Figure  9  on  page  30,  and  Figure  10  on  page  31.  Dashed  lines  indicate  the  true  values 
of  the  parameters.  Solid  lines  show  the  values  the  algorithm  calculated. 

Table  4  on  page  32  contains  the  algorithm's  best  estimates  of  the  parameters, 
along  with  the  number  of  iterations  required  to  converge  to  those  estimates.  It  also 
show's  the  absolute  and  percent  errors  from  the  true  parameters.  For  the  second-order 
test  cases,  convergence  to  the  true  parameter  values  occurred  within  20  iterations.  The 
third-order  test  case  took  14  iterations  to  converge  to  its  steady-state  values;  however, 

these  values  were  not  the  true  parameter  values. 
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Figure  10.  Third-order  test  case  T3.  (A)  MA  parameters.  (B)  AR  parameters. 
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Table  4.  PARAMETER  ESTIMATES  BY  THE  ITERATIVE  MULTICHANNEL 
ALGORITHM. 


TEST 

FILTER 

PA  RA  METER 
ESTIMATE 

ABSOLLTE 

ERROR 

PERCENT 

ERROR 

ITERATIONS 

T2 

0.500 

0.0 

0.0 

5 

-0.393 

-0.007 

1.75 

0.0 

0.0 

!  l.ooo 

0.0 

0.0 

+  0.001 

0.11 

0.247 

-0.003 

1.20 

T2N 

■nnB| 

0.0 

0.0 

20 

+  0.006 

-0.002 

0.0 

0.0 

+  0.004 

0.45 

0.245 

-0.005 

2.00 

T3 

-0.0001 

0.65 

14 

+  0.0025 

5.41 

i  0.0590 

+  0.0128 

27.71 

+  0.0133 

86.36 

0.99 

-0.0 1 

1.0 

-1.79 

+  0.20 

10.05 

1.267 

-0.305 

19.40 

-0.3086 

+  0.1497 

32.66 

2.  Stopping  Parameter 

In  tests  of  third-order  systems,  the  parameter  estimates  swung  through  or  close 
to  the  true  coefficient  values  and  converged  to  values  somewhat  removed  from  the  true 
values.  We  developed  a  stopping  parameter  to  flag  the  iteration  when  the  estimates  were 
closest  to  the  true  values.  This  occurs  when  the  error  is  smallest.  Referring  to 
Figure  7  on  page  25.  if  D(r)  is  equal  to  A(z)  and  C(z)  is  equal  to  B(r),  x  and  z  will  both 
equal  A(z)L'(-)-  The  error  will  be  zero.  The  farther  removed  D(r)  and  C(z)  are  from 
A(r)  and  B(r),  respectively,  the  larger  the  error  becomes. 

The  stopping  parameter  is  calculated  from  the  difference  of  the  z  and  x  values 
at  every’  iteration.  After  the  parameter  vectors  d  and  c  are  estimated  for  a  particular  it¬ 
eration,  the  stopping  parameter  is  calculated  from: 

«*(«)  =  **(«)  -  •*■*(")  =  y kck  -  (4.25) 
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where  k  is  the  iteration  number  and  d  ,  u,  c  and  y  are  defined  in  equations  (4.5)  to  (4.S). 
The  parameter  vectors  d  and  c  represent  the  systems  D(z)  and  C(z),  respectively. 

Figure  11  on  page  34  and  Figure  12  on  page  35  graph  the  stopping  parameter 
(dotted  line)  along  with  the  estimated  and  true  values  of  the  parameters.  They  show  that 
when  the  stopping  parameter  is  smallest,  the  parameters  are  closest  to  their  true  values. 
Table  5  shows  the  improvement  in  the  estimates  of  the  parameters  resulting  from 
choosing  those  values  when  the  stopping  parameter  is  smallest. 


Table  5.  PARAMETER  ESTIMATES  WHEN  STOPPING  PARAMETER  IS 
SMALLEST. 


TEST 

FILTER 

PARAMETER 

ESTIMATE 

ABSOLUTE 

ERROR 

PERCENT 

ERROR 

ITERATIONS 

T3 

■ 

+  0.0002 

1.30 

10 

+  0.0016 

3.46 

0.0536 

+  0.0074 

16.02 

0.0196 

+  0.0042 

27.27 

0.97 

-0.03 

3.0 

-1.S4 

+  0.15 

7.54 

1.375 

-0.197 

12.53 

-0.3672 

+  0.0911 

19.88 

The  stopping  parameter  can  be  used  in  a  real  modeling  situation  because  it 
comes  from  the  data  and  the  estimates  of  the  parameters.  Another  measure  of  how  well 
the  estimates  of  the  parameters  lit  the  actual  system  is  the  norm  of  the  coefficient  error. 
This  cannot  be  used  in  a  real  modeling  situation  however,  because  the  values  of  the  true 
parameters  are  not  known.  We  calculated  it  for  the  test  cases  as  a  check  on  the  appro¬ 
priateness  of  using  the  stopping  parameter.  Figure  13  on  page  36  and  Figure  14  on 
page  37  graph  the  stopping  parameter  (dotted  line)  and  the  norm  of  the  coefficient  error 
for  test  cases  T2  and  T3.  On  both  graphs  the  two  curves  correspond  well.  Both  reach 
their  minimum  value  at  the  same  point,  the  point  where  the  estimates  of  the  parameters 
are  closest  to  their  true  values. 

3.  Linear-prediction  of  the  Denominator  Coefficients 

The  iterative  algorithm  detailed  in  equations  (4.22)  to  (4.24)  starts  by  initializing 
the  AR  parameter  estimates  to: 

c,0)  =  [l  0  ...  0]r  (4.26) 
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Figure  1 1.  Stopping  parameter  example  for  test  case  T2. 
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STOPPING  PARAMETER  STOPPING  PARAMETER 


Figure  13.  Stopping  parameter  and  norm  of  the  coefficient  error  for  test  case  T2. 


where  c0  is  known  to  be  1  from  the  Z  transform  of  the  ARMA  difference  equation, 
equation  (2.1).  The  Z  transform  is  given  by: 


+  c2r  - 

+  -•  ]  =  t(z)[J0  +  d};  1  +  d2z  2  +  •  ] 

(4.27) 

Viz) 

c/q  +  +  d-\JZ  + 

(4-28) 

l\z) 

I  +  C|Z_ 1  +  c2z~2  +  - 

The  initial  estimates  for  the  other  AR  parameters  are  zero.  This  can  be  far  from 
their  actual  values.  A  closer  estimate  of  the  other  AR  parameters  should  result  in 
quicker  convergence  for  all  parameters.  A  closer  estimate  of  the  AR  parameters  can  be 
obtained  by  using  linear-prediction  techniques.  Figure  15  on  page  39  shows  the  ap¬ 
proach  used.  In  Figure  15  on  page  39,  y{n )  is  the  output  from  the  unknown  system. 
The  system  C'(z),  which  is  represented  by  the  vector  c' ,  is  the  linear-prediction  filter  used 
to  estimate  y{n).  It  uses  the  previous  n  —  .V  samples  of  the  output  to  generate  a  current 
estimate.  This  is  given  by  the  equation: 
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Figure  14.  Stopping  parameter  and  norm  of  the  coefficient  error  for  test  case  T3. 

;•(„)  =  yr(«  -  l)c'  (4.29) 

where  c'  is  a  vector  of  the  tap  weights  of  the  autoregressive  process  given  by: 

c'  -  [c'o  c',  ...c'„_A]r  (*4-30) 

and  y(n  —  1)  is  a  vector  of  the  past  output  values  given  by: 

y(«-  I)*  W"-  l)  An~  2)  ...  ><"-‘V)]r  (-4-31) 

Following  least-squares  techniques,  we  form  the  error  between  the  estimate  and  the  ac¬ 
tual  value  of  the  output.  The  sum  of  the  squares  of  the  errors  becomes  the  performance 
criterion.  This  is  differentiated  with  respect  to  the  tap  weights  and  set  equal  to  zero. 
Solving  this  for  the  tap  weights  results  in  the  equation: 

c'-R-’r,,  (4.32) 
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This  is  the  standard  Weiner  filter  solution  [Ref.  8:  p.  32].  It  tells  us  that  the  best  estimate 
of  the  AR  parameters  can  be  found  from  the  correlations  of  the  output  data.  1  ne  matrix 
R>}  is  the  autocorrelation  matrix  of  the  past  outputs  and  is  autocorrelation  vector  of 
the  past  outputs  with  the  current  output.  In  all  cases  tested,  we  did  not  achieve  any 
significant  improvement  in  the  accuracy  of  the  estimates  of  the  parameters,  or  in  the 
speed  of  convergence,  using  the  straight  linear  prediction  of  equation  (4.32). 

A  modification  to  this  approach,  which  we  refer  to  as  modified  linear-prediction, 
uses  correlation  lags  beginning  on  the  order  of  the  MA  portion  of  the  ARM  A  process. 
For  example,  correlations  for  calculating  Rw  for  a  third-order  system  would  start  at  a  lag 
of  three  and  increase  to  a  lag  of  five.  Correlations  for  calculating  rw  would  start  at  a  lag 
of  four  and  increase  to  a  lag  of  six.  This  ensures  that  the  effect  of  the  MA  part  of  the 
unknown  system  is  removed  from  the  linear-prediction  of  the  AR  part.  This  modified 
method  of  linear-prediction  is  given  by  the  equation: 


C'o 

rvy(q)  ryy{q  -  1) 

•  ryy{q-p  +  1) 

ryy{q+  1) 

c'l 

_ 

ryy(q  +  1)  ryy(q) 

•  ryy{q~p  +  2) 

ryy{q  +  2) 

ryy{q  +  p-\)  ryy(q  +  p-  2)  . 

ryy(q) 

ryy{q  +  p) 

where  q  is  the  order  of  the  MA  portion  and  p  is  the  order  of  the  AR  portion.  (Ref.  2: 

p.  182] 

Figure  16  on  page  40  shows  the  results  of  using  modified  linear-prediction  with 
third-order  test  case  T3.  When  comparing  this  graph  to  the  estimates  obtained  without 
linear-prediction,  shown  in  Figure  12  on  page  35,  note  that  the  vertical  axes  have  dif¬ 
ferent  scales.  Table  6  on  page  39  lists  the  values  of  the  estimates  at  iteration  10  and 
compares  them  with  the  true  values  via  the  absolute  and  percent  errors.  A  comparison 
of  Table  6  on  page  39  with  Table  5  on  page  33,  the  best  estimates  without  the  use  of 
modified  linear-prediction,  shows  that  modified  linear  prediction  has  significantly  in¬ 
creased  the  accuracy  of  the  AR  estimates  at  the  tenth  iteration.  The  accuracy  of  the 
MA  estimates  remains  approximately  the  same.  The  tenth  iteration  was  chosen  as  the 
point  to  select  the  parameter  values  because  in  both  cases  this  was  the  iteration  where 
the  stopping  parameter  had  the  smallest  value. 
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Table  6.  PARAMETER  ESTIMATES  USING  MODIFIED 

LINEAR-PREDICTION  FOR  INITIAL  ESTIMATE  OF  AR  PARAME¬ 
TERS. 


TEST 

FILTER 

PARAMETER 

ESTIMATE 

ABSOLUTE 

ERROR 

PERCENT 

ERROR 

ITERATIONS 

T3 

+  0.0002 

10 

. .  ;  •  .  ■ 

+  0.0023 

+  0.0050 

10.S0 

0.0101 

+  0.0053 

34.42 

1.00 

0.0 

0.0 

-1.98 

+  0.01 

0.50 

1.553 

-0.019 

1.21 

-0.4458 

+  0.0125 

2.73 

D.  FORMULATION  OF  THE  SEQUENTIAL  MULTICHANNEL  APPROACH 

To  decrease  the  computational  intensity  of  updating  the  estimates  of  the  AR  and 
MA  parameters  due  to  new  data,  an  algorithn  to  sequentially  process  the  data  has  been 


Figure  16.  Parameter  estimates  for  test  case  T3  using  modified  linear-prediction  for 
initial  estimate  of  AR  parameters. 
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developed.  Development  begins  with  the  performance  criterion  seen  previously  for  the 
block  data  case  in  equation  (4.10).  From  that  starting  point,  equations  are  developed 
which  relate  new  estimates  of  the  parameters  to  the  previous  estimates  and  the  new  data. 
Separate  but  similar  equations  are  developed  for  the  \1A  and  the  AR  coefficients.  Due 
to  the  similar  nature  of  the  development  of  these  equations,  only  the  development  of  the 
equations  for  the  AR  coefficients  is  presented  here.  However,  the  final  results  for  both 
the  AR  and  MA  coefficients  are  given. 

The  performance  criterion  can  be  written  as: 


n  n 

J  -  “ 

<=0  /=0  /= 0 


(4.34) 


Expanding  this  results  in: 


J 


n 

-  2z^c  +  c7W« 

i=0 


(4.35) 


where  c  and  y  are  defined  in  equations  (4.7)  and  (4.8)  and  z  is  the  scalar  signal  at  the 
output  of  D.  Differentiating  the  performance  criterion  with  respect  to  c  and  setting  the 
result  equal  to  zero  yields: 


cJ_ 

cc 


=  0  = 


Equation  (4.36)  can  also  be  written  as 


Z*,-(2V)e 

1=0  1=0 

Solving  for  the  AR  parameter  vector  results  in: 

« = 

1=0  1=0 


(4.36) 


(4.37) 


(4.38) 
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Because  c  is  an  estimate  based  on  data  available  through  time  n  we  signify  this  by  in¬ 
troducing  the  index  n  on  c  to  yield: 


n  n 

/=o  1=0 

The  first  step  in  formulating  the  sequential  algorithm  is  to  develop  an  update 
equation  for  the  estimate  of  the  AR  parameters.  Applying  the  method  presented  in 
Graupe  [Ref.  9:  p.  124J,  we  first  define  a  new  matrix  P;'  as 

n 

p;'  =  Xy,y,T  (4.40) 

1=0 

This  is  a  matrix  of  the  output  data  of  the  unknown  system.  At  the  previous  time  n  —  1 
this  matrix  is  written  as: 


n- 1 

p;i, =£>•/>/'  <‘>41) 

(=0 

By  substituting  equation  (4.40)  into  equation  (4.39)  the  estimate  of  the  AR  parameters 
can  be  rewritten  as: 


C*  =  p (4.42) 

i= 0 

The  right  side  of  equation  (4.42)  needs  to  be  converted  into  an  expression  containing  the 
previous  estimate  of  the  parameters  plus  a  correction  term.  It  needs  the  past  value  of 
c„  which  is  c„_,.  From  equation  (4.39)  we  can  write: 


c 


n— 1 


(Iy,yrr'I>< 


1=0  1=0 


(4-43) 


This  can  be  rewritten  as: 


n-l  n—1 


i=0  /« 0 


(4.44) 
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Premultiplying  equation  (4.42)  by  P;1  and  then  separating  the  last  term  from  the  rum¬ 
ination  results  in: 


n— 1 

p« ’c« = Yjz‘y‘ + 2«>* 

z=0 


(4-45) 


Substituting  for  Xzfy,  in  equation  (4.45)  from  its  equivalent  expression  in  equation  (4.44) 

i=0 

yields: 


n—  i 

lc  *  1  +  (4.46) 

i=0 

By  adding  yBy^c„..,  to  the  right  side  of  equation  (4.46)  and  grouping  it  with  the  summa¬ 
tion,  the  summation  will  range  from  /  =  0  to  n  .  In  order  not  to  affect  the  value  on  the 
right  side  of  equation  (4.46),  y„y„rcn_,  must  also  be  subtracted  from  the  right-hand  side, 
which  yields: 


n-1 

p;i,  c  =  (^y(y/)c„_,  +  zn  y„  +  y„yJcB_l  -  yny[c„_,  (4.47) 

i*0 

Combining  y„yfo_,  with  the  summation  as  describe  above  results  in: 

n 

Pn-lC  =  +  Zn}'n  ~  y^Un-l  (4-48) 

i=0 

rt 

Replacing  £y,y f  with  its  equivalent  expression  from  equation  (4.40)  yields: 

i=l' 

-  p;1'.-,  +  yJf.  -  K49) 

Premultiplying  by  P„  results  in  the  following  equation  for  the  update  of  the  estimate  of 
the  AR  parameters: 

c*  =  c„_,  +  Pnyn(z„  -  y^n-i)  (4-5°) 

This  is  the  desired  result.  It  relates  the  estimate  of  the  parameters  at  time  Ar  to  the  es¬ 
timate  at  the  previous  time,  AT  —  1,  plus  the  new  output  data  vector,  y„,  and  the  error  at 
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time  N.  Error  is  represented  by  the  z„  —  yfo.,  term.  The  corresponding  equation  for  the 
MA  parameters  is: 

d»  -  d„_,  +  Qnun(x„  -  ujd„_,)  (4.51) 

In  equation  (4.51),  Q  is  a  matrix  of  the  input  data  of  the  unknown  system  given  by: 

n 

Q»"'  -  X"Ar  (4.52) 

i=Q 

Finally,  we  need  a  sequential  update  for  P„  and  Q„  to  complete  the  sequential  algo¬ 
rithm.  This  is  accomplished  by  using  a  form  of  the  matrix  inversion  lemma. 

By  pulling  the  last  term  out  of  the  summation,  the  definition  of  P;1  given  in  equation 
(4.40)  can  be  rewritten  as: 

p;'=2y'y'r+y^r 

/= 0 

n- 1 

Substituting  for  £y,y J  its  equivalent  expression  from  equation  (4.41)  results  in 

f=0 

p;'  -  p;',  +  y,y,r 

Inverting  both  sides  of  the  equation  results  in: 

p„  =  (P„li  +  ynyl)~'  (4.55) 

Let  A  =  P;i„  B  =  y„  C  =  1  ,  and  D  =  yj.  Then,  by  the  matrix  inversion  lemma: 

P„  -  A-1  -  A_1B(C'’  +  DA-IB)-1DA-1  (4.56) 

Substituting  the  appropriate  expressions  into  equation  (4.56)  results  in: 

p„ = (p;' ,r'  -  (p;i,r'j»[i + ySp;i,rly»r1yJ'(p;l,r'  (4.57) 

This  reduces  to: 

p*  =  p„_,  -  P„_,y„[i  +  yX-.yJ-'y^-,  (4-58) 

Using  this  same  procedure,  the  update  for  Q„  is: 


(4.53) 


(4.54) 
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(4.59) 


Q*  =  Q/,-1  -  Q«-1««[1  +  u^iUj-’u^Q^, 

A  reduction  in  the  computational  intensity  has  been  achieved  by  reducing  the  matrix 
inversion  of  equation  (4.55)  to  inversion  of  a  scalar  in  equation  (4.57).  Inversion  of 
these  scalars  is  much  simpler  than  inversion  of  the  and  R.,,  matrices  of  the  block 
processing  case. 

The  sequential  multichannel  algorithm  is  summarized  below: 


The  parameter  update  equations: 

c*  =  c„_,  +  P„y „{zn  -  y[cn_,) 

(4.60) 

=  d„-i  +  Q  *»*„(*„  -  uX-i) 

(4.61) 

The  update  equations  for  the  P  and  Q  matrices: 

P*  =  P„_,  -  P„-,y„[i  +  y„rPrt_,y„]_1ynrP„_, 

(4.62) 

Qn  =  Q„-1  -  Q*-iU„[I  +  u„rQn.1uJ'’unrQ„_1 

(4.63) 

The  reduction  in  computational  intensity  comes  with  a  trade-off.  Now  the  algo¬ 
rithm  is  more  complex  to  use.  Updates  are  required  for  P  and  Q  as  well  as  c  and  d  where 
before,  in  the  block  multichannel  algorithm,  updates  were  only  required  for  c  and  d.  But, 
as  in  the  sequential  IV  case,  an  added  advantage  of  the  sequential  multichannel  algo¬ 
rithm  is  it  allows  updates  of  the  estimates  of  the  parameters  based  upon  new  data. 
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V.  SUMMARY 


In  this  thesis  we  set  out  to  develop  two  algorithms  for  modeling  unknown  systems 
as  ARMA  processes.  These  are  the  IV  method  of  system  identification  presented  in 
Chapter  3,  which  is  a  modification  of  the  method  of  ordinary  least-squares,  and  the  it¬ 
erative  multichannel  method  presented  in  Chapter  4. 

The  IV  method  is  not  a  new  concept  in  either  its  block  or  sequential  processing 
forms.  However,  our  derivation  of  the  sequential  algorithm  was  done  independently  of 
other  IV  sequential  algorithms.  We  chose  the  IV  method  because  it  reportedly  has  good 
noise  handling  capabilities  and  yields  consistent  and  unbiased  estimates  of  the  unknown 
system  s  parameters.  It  also  remains  as  easy  to  use  as  the  method  of  ordinary'  least- 
squares.  We  also  wanted  to  gain  familiarity  with  it  because  it  Was  a  possible  candidate 
for  incorporation  into  the  multichannel  method. 

Operating  alone,  the  IV  method  produces  accurate  estimates  of  the  unknown  sys¬ 
tem's  parameters.  Convergence  was  within  20  iterations  for  several  second-order  sys¬ 
tems  that  we  tested.  Convergence  slows  down  as  the  system  order  increases.  However, 
the  results  do  converge  to  the  actual  system  parameters  given  sufficient  processing  time. 
The  performance  of  the  IV  method  is  similar  to  the  performance  of  the  method  of  ordi¬ 
nary  least-squares. 

The  proposed  iterative  multichannel  algorithm  is  new  in  both  its  block  and  sequen¬ 
tial  processing  forms  to  the  best  of  our  knowledge.  It  is  very  simple  to  use  in  the  block 
form.  It  achieves  accurate  results  for  second-order  systems  but  worse  results  for  third- 
order  systems  with  block  correlation  elements  calculated  based  on  only  500  data  points. 
Implementing  the  stopping  parameter  increases  the  accuracy  when  the  parameter  esti¬ 
mates  converge  but  not  to  the  true  parameters.  Due  to  its  ability  to  separately  process 
the  input  and  output  data  from  the  unknown  system,  correlation  matrices  in  the  multi¬ 
channel  block  processing  case  are  half  the  size  of  correlation  matrices  required  for  the 
single  channel  block  processing  case.  This  feature  reduces  the  computational  intensity 
over  what  is  required  for  the  conventional  least-squares  processing  case.  The  number 
of  iterations  required  for  convergence  seems  to  be  independent  of  the  order  of  the  sys¬ 
tem.  However,  the  accuracy  of  the  estimates  suffer  as  the  order  of  the  system  increases. 
Using  linear  prediction  to  estimate  the  initial  values  of  the  AR  parameters  did  not  speed 
up  convergence  or  increase  the  accuracy  of  the  parameter  estimates.  However,  using 
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modified  linear  prediction  significantly  increased  the  accuracy  of  the  AR  parameter  es¬ 
timates,  although  it  had  no  effect  on  the  MA  parameter  estimates. 

We  formulated  the  multichannel  sequential  algorithm.  This  allows  the  estimates  of 
the  parameters  to  be  updated  as  new  data  becomes  available.  But  we  have  not  tested 
this  algorithm.  It  needs  checking  using  a  variety  of  second  and  third-order  test  cases  to 
verify  that  it  works.  During  testing,  guidelines  need  to  be  developed  for  the  best  meth¬ 
ods  to  initialize  the  P  and  Q  matrices  to  achieve  the  quickest  convergence  and  the  most 
accurate  parameter  estimates. 

As  mentioned  above,  one  of  the  reasons  for  investigating  the  IV  method  of  system 
identification  was  for  possible  inclusion  into  the  multichannel  algorithm,  the  hope  being 
that  the  favorable  noise  performance  of  the  IV  method  would  improve  the  performance 
of  the  multichannel  method.  This  is  another  area  that  remains  unexplored. 

The  block  multichannel -and  IV  methods  achieved  similar  results  for  second-order 
test  cases.  Convergence  to  the  actual  system  parameters  came  within  20  iterations  for 
both  algorithms.  However,  for  third-order  systems,  convergence  was  much  quicker  with 
the  multichannel  block  algorithm  than  with  the  sequential  IV  algorithm.  But  the  pa¬ 
rameter  estimates  by  the  IV  method  were  more  accurate  than  by  the  multichannel  block 
method.  A  combination  of  the  two  algorithms  has  the  potential  for  incorporating  their 
unique  advantages  into  a  better  overall  parameter  estimation  method. 

Areas  for  further  research  are  listed  below: 

•  Verify  that  the  multichannel  sequential  algorithm  developed  in  Chapter  4  works  as 
a  means  of  modeling  an  unknown  system  as  an  ARM  A  process. 

•  Investigate  the  possibility  of  incorporating  the  IV  method  into  the  multichannel 
sequential  algorithm. 

•  Analyze  why  initializing  the  AR  parameters  to  the  values  calculated  by  linear  pre¬ 
diction  improves  the  speed  of  convergence  of  the  AR  parameters  in  the  multi¬ 
channel  block  algorithm  but  does  not  improve  the  convergence  of  the  MA 
parameters.  Identify  a  method  for  obtaining  an  initial  estimate  of  the  MA  param¬ 
eters  to  improve  their  speed  of  convergence  and  accuracy. 

•  Investigate  the  effects  of  increasing  the  number  of  data  points  used  to  calculate  the 
correlation  matrices  for  the  multichannel  block  algorithm  on  the  accuracy  of  the 
parameter  estimates  and  their  speed  of  convergence. 

•  Investigate  the  performance  of  the  IV  method  with  noise  present.  Compare  this 
performance  with  the  performance  of  the  method  of  ordinary  least-squares  with 
noise  present. 
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APPENDIX 
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A.  INSTRUMENTAL  VARIABLE  ALGORITHM 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


PAUL  DAL  SANTO 
IV  ALGORITHM 
4/12/88 

THIS  PROGRAM  CALCULATES  THE  AR  AND  MA  PARAMETERS  OF  A 
TEST  SYSTEM  BASED  UPON  ITS  INPUT  AND  OUTPUT  DATA 
BY  USING  THE  SEQUENTIAL  IV  METHOD. 


******* 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


**** 


VARIABLE  DEFINITIONS 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


* 

* 

* 

* 

* 

* 

* 


RIVCOF  ARRAY  FOR  STORING  THE  AR  AND  MA  PARAMETERS 
THE  PROGRAM  CALCULATES 

Z  VECTOR  OF  DATA  FROM  THE  OUTPUT  OF  THE  AUXILIARY 

MODEL  AND  THE  INPUT  TO  THE  TEST  SYSTEM 
Z  TPO  TRANSPOSE  OF  VECTOR  Z 

X  VECTOR  OF  DATA  FROM  THE  OUTPUT  AND  INPUT  OF  THE 

TEST  SYSTEM 

X  TPO  TRANSPOSE  OF  THE  X  VECTOR 

U  STORAGE  FOR  INPUT  DATA 

Y  STORAGE  FOR  OUTPUT  OF  TEST  SYSTEM 

W  STORAGE  FOR  OUTPUT  OF  THE  AUXILIARY  MODEL 

QMAT  THE  Q  MATRIX  OF  THE  IV  ALGORITHM 

IV  VECTOR  OF  PARAMETERS  CALCULATED  BY  THE  ALGORITHM 

COR  RESULT  OF  INTERMEDIATE  STEP  IN  ALGORITHM  CALCULATION 

COEF  VECTOR  OF  TRUE  PARAMETERS  OF  TEST  SYSTEM 

SCALAR  RESULT  OF  SCALAR  INVERSION  IN  INTERMEDIATE  STEP  COR 
SCALR2  INTERMEDIATE  STEP  WHEN  CALCULATING  THE  NEW  IV  VECTOR 
SEED  INITIALIZATION  FOR  IMSL  GAUSSIAN  ROUTINE 

WSIZE  ORDER  OF  AR  PART  OF  THE  AUXILIARY  MODEL 

YSIZE  ORDER  OF  THE  AR  PART  OF  THE  TEST  SYSTEM 

USIZE  ORDER  OF  THE  MA  PART  OF  THE  TEST  SYSTEM  AND  AUXILIARY 

MODEL 


VARIABLES  THAT  END  IN  R  CONTAIN  THE  ROW  SIZE  OF  THEIR  RESPECTIVE 
MATRICES.  VARIABLES  THAT  END  IN  C  CONTAIN  THE  COLUMN  SIZE  OF 
THEIR  RESPECTIVE  MATRICES. 

VARIABLE  DECLARATIONS  **** 

COMMON  /D/  RIVCOFCO: 1000,10} 

REAL  Z(10,10),Z  TPO(  10 , 10)  ,X(  10, 10)  ,X  TPO(IO.IO') 

REAL  U( 1) ,Y( 10, 10) ,W( 10,10) 


IVA00010 

IVA00020 

IVA00030 

IVA00040 

IVA00050 

IVA00060 

IVA00070 

IVA00080 

IVA00090 

IVA00100 

IVA00110 

IVA00120 

IVA00130 

IVA00140 

IVA00150 

IVA00160 

IVA00170 

IVA00180 

IVA00190 

IVA00200 

IVA00210 

IVA00220 

IVA00230 

IVA00240 

IVA00250 

IVA00260 

IVA00270 

IVA00280 

IVA00290 

IVA00300 

IVA00310 

IVA00320 

IVA00330 

IVA00340 

IVA00350 

IVA00360 

IVA00370 

IVA00380 

IVA00390 

IVA00400 

IVA00410 

IVA00420 

IVA00430 

IVA00440 

IVA00450 

IVA00460 

IVA00470 

IVA00480 


REAL  QMAT( 10 , 10)  ,Q  TEMP( 10 , 10) ,QTEMP2( 10 , 10) ,QTEMP3( 10, 10) 
REAL  IV( 10,10), IVTEMP( 10,10), COR( 10,10), COEF( 10,10) 

REAL  SCALAR, SCALR2 

DOUBLE  PRECISION  SEED 

INTEGER  I , J,K,WSIZE , YSIZE ,USIZE 

INTEGER  ZMR,ZMC,ZTR,ZTC,XMR,XMC,XTR,XTC 

INTEGER  UMR , UMC , YMR , YMC , WMR , WMC 

INTEGER  QMR , QMC , QTR , QTC , QT2R , QT2C , QT3R , QT3C 

INTEGER  I VR , I VC , I VTR , I VTC , CMR , CMC , COEFMR , COEFMC 

INTEGER  IVCR, IVCC 

LOGICAL  EOF 

READ(4,*,END=22)  YSIZE, US I ZE, COEFMR, COEFMC, ITERA 
CALL  RDMAT( COEF , COEFMR , COEFMC ) 

*  INITIALIZE  VARIABLES 

EOF  =  . FALSE. 

XTR  =  COEFMC 
XTC  =  COEFMR 
ZTR  =  COEFMC 
ZTC  «  COEFMR 
IVR  =  COEFMR 
IVC  =  COEFMC 
QMR  =  COEFMR 
QMC  =  COEFMR 
SEED  =  888. 0D1 
IVCR  =  0 
IVCC  =  COEFMR 
WSIZE  ■  YSIZE 

*  ZERO  OUT  THE  IV  PARAMETER  VECTOR,  THE  AUXILIARY  MODEL  DATA  VECTOR 

*  AND  THE  TEST  SYSTEM  DATA  VECTOR. 

CALL  INIT(IV,IVR,IVC,0. 0) 

CALL  INIT(Z  TP0,ZTR,ZTC,0.  0) 

CALL  INIT(X  TPO, XTR, XTC, 0.0) 

*  INITIALIZE  THE  QMAT  AS  A  DIAGONAL  MATRIX  WHOSE  DIAGONAL  ELEMENTS 

*  EQUAL  100 

CALL  INITDCQMAT.QMR, QMC, 100.0) 

*  GET  VALUE  FOR  U(0).  U  IS  A  GAUSSIAN  RANDOM  VARIABLE. 

CALL  GGNML  (SEED,1,U) 

*  SHIFT  U(0)  INTO  X  &  Z  VECTORS  TO  CREATE  X(0)  &  Z(0) 

Y(  1 , 1)  =0.0 

CALL  SHIFT(X  TPO, XTC .YSIZE ,USIZE,Y( 1 , 1) ,U( 1)) 

W(l,l)  =  0.0 

CALL  SHIFT(Z  TPO, ZTC , WSIZE ,USIZE,W( 1 , 1) ,U( 1) ) 


1VA00490 

IVA00500 

IVAOOSiO 

IVA00520 

IVA00530 

IVA00540 

IVA00550 

IVA00560 

IVA00570 

IVA00580 

IVA00590 

IVA00600 

IVA00610 

IVA00670 

IVA00680 

IVA00690 

IVA00700 

IVA00630 

IVA00640 

IVA00650 

IVA00710 

IVA00720 

IVA00730 

IVA00740 

IVA00750 

IVA00760 

IVA00770 

IVA00780 

IVA00790 

IVA00800 

IVA00810 

IVA00820 

IVA00830 

IVA00840 

IVA00850 

IVA00860 

IVA00870 

IVA00880 

IVA00890 

IVA00900 

IVA00910 

IVA00920 

IVA00930 

IVA00940 

IVA00950 

IVA00960 

IVA00970 

IVA00990 

IVA01000 

IVA01010 

IVA01020 

IVA01030 

IVA01040 

IVA01050 
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*  CALCULATE  Y(0)  =  X  TPO(O)  *  COEFFICIENT  VECTOR 

CALL  MULTI (X  TPO,XTR,XTC,COEF,COEFMR,COEFMC,Y, 1 , 1) 

*  CALCULATE  W(O)  =  2  TPO(O)  *  IV  VECTOR 

CALL  MULTI (Z  TPO,ZTR,ZTC,IV,IVR,IVC,W,l,l) 

DO  90  I  =  0 , ITERA 

CALL  GGNML(SEED, 1,U) 

CALL  TPOSE( IV , IVR , IVC , IVTEMP , IVTR , I VTC) 

*  SAVE  THE  IV  PARAMETERS 

DO  91  J  =  l.IVTC 

RIVCOF( I ,  J)  «  IVTEMPC 1 ,  J) 

91  CONTINUE 

CALL  PRMAT( IVTEMP , IVTR , IVTC ) 

*  SHIFT  Y(M)  AND  U(M+1)  INTO  X  TPO(M)  TO  GET  X  TP0(M+1) 

CALL  SHIFT(X  TPO.XTC , YSIZE .USIZE ,Y( 1 , 1) ,U( 1)) 

*  SHIFT  W(M)  AND  U(M+1)  INTO  Z  TPO(M)  TO  GET  Z  TP0(M+1) 

CALL  SHIFT(Z  TPO.ZTC ,WSIZE, USIZE, W( 1 , 1) ,U( 1)) 

*  CALCULATE  Y(M+1)  AND  Z(M+1) 

CALL  MULTI (X  TP0,XTR,XTC,C0EF,C0EFMR,C0EFMC,Y,1,1) 

CALL  MULTI (Z  TPO,ZTR,ZTC , IV, IVR, IVC,W, 1 , 1) 

*  CALCULATE  THE  NEW  Q  MATRIX 

CALL  MULTI (X  TPO,XTR,XTC,QMAT,QMR,QMC,Q  TEMP,QTR,QTC) 

CALL  TPOSE( Z  TPO,ZTR,ZTC,Z,ZMR,ZMC) 

CALL  CORE( ('MAT , QMR , QMC , Z , ZMR , ZMC , X  TPO,XTR,XTC,COR,CMR,CMC) 
CALL  MULTI ( COR, CMR, CMC, Q  TEMP , QTR , QTC , QTEMP2 , QT2R , QT2C ) 

CALL  SUBTRC ( QMAT , QMR , QMC , QTEMP2 , QT2R , QT2C , QMAT , QMR , QMC ) 

*  CALCULATE  THE  NEW  IV  VECTOR 

CALL  MULTI (X  TPO,XTR,XTC, IV, IVR, I VC, IVTEMP, IVTR, IVTC) 

SCALR2  =  Y{ 1,1)  -  IVTEMPC 1 , 1) 

CALL  SMULTI ( SCALR2 , COR , CMR , CMC ) 

CALL  ADD( I V , IVR , IVC , COR , CMR , CMC , IV , IVR , IVC ) 

90  CONTINUE 

*  PLOT  THE  IV  PARAMETERS  VS  THE  ITERATION  NUMBER 

CALL  PL0T2( ITERA,USIZE,YSIZE ,COEF) 

22  STOP 

END 

*  *****************  ** ** ****** *** ********************* ************* 
SUBROUTINE  CORE ( MAT 1 , I 1R , I 1 C , MAT2 , I2R, I2C, MAT3 , I 3R , I 3C , RMAT , IRR , 

+IRC) 


IVA01070 
IVA01080 
IVA01090 
IVA01100 
IVA01110 
IVA01 120 
IVA01130 
IVA01140 
IVA01160 
IVA01180 
IVA01200 
IVA01210 
IVA01220 
IVA01230 
IVA01240 
IVA01250 
1VA01260 
IVA01280 
IVA01290 
IVA01300 
IVA01320 
IVA01330 
IVA01340 
IVA01360 
IVA01370 
IVA01380 
IVA01400 
IVA01410 
IVA01420 
IVA01430 
IVA01440 
IVA01450 
IVA01460 
IVA01470 
IVA01480 
IVA01490 
IVA01500 
IVA01510 
IVA01520 
IVA01530 
IVA01540 
IVA01550 
IVA01560 
IVA01570 
IVA01590 
IVA01600 
IVA01610 
IVA01630 
IVA01640 
IVA01650 
IVA01660 
IVA01670 
IVA01930 
IVA01940 
IVA01950 
IVA01960 


REAL  MAT1( 10,10) ,MAT2( 10,10) ,MAT3( 10,10) ,RMAT( 10,10) 

REAL  Q  TEMP( 10,10) ,QTEMP2( 10,10) 

INTEGER  IRR , IRC , QTR , QTC , QT2R , QT2C 

*  CALCULATE  THE  CORE:  Q(M)Z(M+1) #1+X' (M+1)Q(M)Z(M+1)  **-l 

*  MAT1  IS  THE  Q  MATRIX,  MAT2  IS  THE  Z  VECTOR,  AND  MAT3  IS  THE 

*  X  VECTOR. 

CALL  MULTI(MAT1,I1R,I1C,MAT2,I2R,I2C,Q  TEMP, QTR, QTC) 

CALL  MULTI(MAT3, I3R, I3C ,Q  TEMP, QTR, QTC, QTEMP2,QT2R,QT2C) 
SCALAR  =  1/(1  +  QTEMP2( 1,1)) 

CALL  SMULTIC SCALAR, Q  TEMP , QTR , QTC ) 

CALL  ADD(Q  TEMP , QTR , QTC , Q  TEMP, QTR, QTC, RMAT, IRR, IRC) 

CALL  SMULTICO.  5 , RMAT, IRR, IRC) 

RETURN 

END 

•> V  •it'kicirirlKjKirlt-k’kft'kicir'kjcirlrk'k'krk'kirkic'klc'k 

SUBROUTINE  PL0T2( ITERA , ICURVN , ICURVD , MAT1 ) 

*  ****************************** 


COMMON  /D/  RIVCOF(0: 1000,10) 

REAL  X( 0:  1000), Y(0:  1000) ,MAT1( 10, 10) .MAX, MIN 
INTEGER  I, J, ITERA, ICURVN, ICURVD, ISTP 

CALL  LI M ITS ( I CURVN , I CURVD , NMAX , NMIN , NSTP , 
+DMAX , DMIN , DSTP , ITERA ) 

*  GENERATE  THE  ITERATION  NUMBER 

DO  90  I  =  0, ITERA 
X(I)  =  I 
Y(  I)  =  0.0 
90  CONTINUE 

*  CALCULATE  X  AXIS  LABELING  INTERVAL 

ISTP  =  ITERA/ 10 

*  SET  UP  THE  PLOT 

CALL  SHERPAC ' IVGRAF  ' , ' A ' , 3 ) 

CALL  RESET('ALL') 

CALL  PAGE(8.511.0) 

CALL  HEIGHT(0. 14) 

CALL  HWROT('AUTO') 

CALL  XINTAX 
CALL  YAXANG(O) 

CALL  PHYSORQ.  5,6.  0) 

CALL  AREA2D(5. 0,3.5) 

CALL  COMPLX 

CALL  XNAME( ' ITERATIONS?' ,100) 

CALL  YNAMEC COEFFICIENT  VALUES?' ,100) 

CALL  MESSAGC(A)?' ,100,2.  4,-0.  8) 

CALL  THKFRM(0.  03) 

CALL  FRAME 
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CALL  GR AF ( 0 , I STP , I TERA , NM I N , N STP , NMAX ) 

*  PLOT  THE  NUMERATOR  VALUES 

DO  93  J  =  ICURVD  +  1,ICURVN  +  ICURVD 
DO  94  I  =  0 , ITERA 
Y( I)  =  RIVCOF( I ,  J) 

94  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

93  CONTINUE 

*  PLOT  DASHED  LINES  FOR  THE  COEFS'  TRUE  VALUES 

CALL  DASH 

*  PLOT  NUMERATOR  COEFS'  TRUE  VALUES 

DO  95  K  =  ICURVD  +  1, ICURVD  +  ICURVN 
DO  96  J  =  0, ITERA 
Y( J)  =  MAT1(K,1) 

96  CONTINUE 

CALL  CURVE ( X , Y , ITERA , 0 ) 

95  CONTINUE 
CALL  ENDGR(O) 

*  SET  UP  SECOND  PLOT  FOR  DENOMINATOR  PARAMETERS 

CALL  RESET(’DASH') 

CALL  PHYSORC  1.5,1.  5) 

CALL  AREA2D(5.  0,3.5) 

CALL  COMPLX 

CALL  XNAME( ' ITERATIONS* ’ , 100) 

CALL  YNAME( ' PARAMETER  VALUES* ' , 100) 

CALL  MESSAG( ' (B)$ ' , 100,2.  4,-0.  8) 

CALL  THKFRM(0.  03) 

CALL  FRAME 

CALL  GRAF( 0 , ISTP , ITERA, DMIN,DSTP,DMAX) 

*  PLOT  THE  DENOMINATOR  VALUES 

DO  91  J  =  1, ICURVD 
DO  92  I  =  0, ITERA 

Y( I)  =  -RIVCOF( I , J) 

92  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

91  CONTINUE 

*  PLOT  DENOMINATOR  COEFS'  TRUE  VALUES 

CALL  DASH 

DO  99  K  =  1, ICURVD 

DO  100  J  =  0, ITERA 
Y( J)  =  -MAT1(K, 1) 

100  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

99  CONTINUE 
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IVA04420 
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IVA04570 

IVA04580 
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IVA04610 

IVA04620 

IVA04630 

IVA04640 
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IVA04770 
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IVA04790 
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on  *  *  *  *00000 


CALL  ENDPL(O) 
CALL  DONEPL 
RETURN 
END 


* 


* 


SUBROUTINE  SUBTRC( MAT1 , I 1R , I 1C , MAT2 ,  I2R,I2C, RMAT , IRR , IRC ) 


****  PURPOSE  -  ROUTINE  SUBTRACTS  MAT2  FROM  MAT1  AND  PUTS  THE 
*  RESULT  IN  A  THIRD  MATRIX. 

REAL  MAT1( 10 , 10) ,MAT2( 10,10) ,RMAT( 10,10) 

INTEGER  IRR, IRC 

CALL  SMULTIC-1.  0 ,MAT2 , 12R, I2C) 

CALL  ADD(MAT1,I1R,I1C,MAT2,I2R,I2C, RMAT, IRR, IRC) 

IRR  =  HR 
IRC  =  I1C 
RETURN 
END 


* 

•k 


SUBROUTINE  TPOSE( MAT 1 , I 1R , I 1C , RMAT , IRR , IRC ) 


****  PURPOSE -SUBROUTINE  TRANSPOSES  A  MATRIX  AND  PUTS  THE  RESULT 

INTO  A  NEW  MATRIX 

REAL  MAT1( 10, 10) ,RMAT( 10 , 10) 

INTEGER  I, J, IRR, IRC 

DO  93  1=1, I1C 
DO  94  J=1 , I1R 

RMAT(I,J)  =  MAT1(J,I) 

94  CONTINUE 

93  CONTINUE 
IRR  =  I1C 
IRC  =  HR 
RETURN 
END 
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B.  MULTICHANNEL  ALGORITHM 


* 


* 


*  PAUL  DAL  SANTO  8/15/88  * 

*  * 


*  TWO-CHANNEL  SYSTEM  IDENTIFICATION  ALGORITHM 


*  * 

*  PROGRAM  CALCULATES  THE  AR  AND  MA  PARAMETERS  * 

*  BASED  UPON  THE  INPUT  AND  OUTPUT  DATA  OF  A  * 

*  TEST  SYSTEM.  * 

*  * 

*  SHIFT  USED  TO  CALC  RYY  FOR  THE  LINEAR  * 


DUA00010 

DUA00020 

DUA00030 

DUA00040 

DUA00050 

DUA00060 

DUA00070 

DUA00080 

DUA00090 

DUA00100 

DUA00110 


nooooon  *  a 


*  OF  THE  AR  PARAMETERS  IS  READ  FROM  * 

*  THE  COEFF  DATA  FILE  * 

*  * 

*  NUMBER  OF  ITERATIONS  IS  READ  FROM  COEFF  * 

*  DATA  FILE  * 

*  * 

*  * 

*  * 


************** k**kk*k*kkk ************************************ 


* 


**** 


VARIABLE  DEFINITIONS 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


* 

* 

* 

* 

* 

* 

* 


RAWDAT  ARRAY  WHICH  CONTAINS  THE  INPUT  AND  OUTPUT  DATA 
OF  THE  SYSTEM  UNDER  TEST 

El  ARRAY  FOR  STORING  THE  STOPPING  PARAMETER  AND  THE 

COEFFICIENT  ERROR 

ARRD  ARRAY  FOR  STORING  THE  AR  PARAMETER  ESTIMATES 

ARRN  ARRAY  FOR  STORING  THE  MA  PARAMETER  ESTIMATES 

RYYM  AUTOCORRELATION  MATRIX  OF  THE  OUTPUT  DATA 
RYUM , RUYM  CROSSCORRELATION  MATRICES 

RUUM  AUTOCORRELATION  MATRIX  OR  THE  INPUT  DATA 
RUUINV  INVERSE  OF  THE  AUTOCORRELATION  MATRIX  OF  THE  INPUT  DATA 

RYYINV  INVERSE  OF  THE  AUTOCORRELATION  MATRIX  OF  THE  OUTPT  DATA 

DM  VECTOR  OF  CURRENT  MA  PARAMETER  ESTIMATES 

CM  VECTOR  OF  CURRENT  AR  PARAMETER  ESTIMATES 

TRUNRM  FUNCTION  WHICH  CALCULATES  THE  NORM  OF  THE  TRUE  VALUES 
OF  THE  TEST  SYSTEM'S  PARAMETERS 
WKMAT  WORKING  MATRIX  FOR  DOING  MATRIX  INVERSIONS 
X  TPO  TRANSPOSE  OF  THE  VECTOR  OF  INPUT  DATA 
Y  MATRIX  FOR  THE  OUTPUT  OF  THE  TEST  SYSTEM 

U  INPUT  TO  THE  TEST  SYSTEM 

COEFM  VECTOR  OF  TRUE  COEFFICIENTS  OF  THE  TEST  SYSTEM 
ITERA  NUMBER  OF  ITERATIONS  TO  PERFORM 

CORRLN  LENGTH  OF  CORRELATIONS  TO  USE  TO  CALCULATE  RYY,  RUU 
RYU,  AND  RUY 

YSIZE  ORDER  OF  THE  AR  PART  OF  THE  TEST  SYSTEM 

USIZE  ORDER  OF  THE  MA  PART  OF  THE  TEST  SYSTEM 

KTIME  THE  CURRENT  ITERATION 

SHFT  SIZE  OF  THE  STARTING  LAG  FOR  LINEAR  PREDICTION  OF 

THE  AR  PARAMETERS 

SEED  INITIALIZATION  PARAMETER  FOR  IMSL  ROUTINE  WHICH 

GENERATES  RANDOM  GAUSSIAN  NUMBERS 


*  INTEGER  VARIABLES  THAT  END  WITH  R  CONTAIN  THE  ROW  SIZE  OF 

*  A  PARTICULAR  ARRAY.  INTEGER  VARIABLES  THAT  END  WITH  C  CONTAIN 

*  THE  COLUMN  SIZE  OF  A  PARTICULAR  ARRAY. 


* 


VARIABLE  DECLARATIONS 


COMMON  /D/  RAWDAT(2,0:  2000), El(2,0:  1000), 

+ARRD(0:  1000,5) ,ARRN(0:  1000,5) 

REAL  RYYM(5 ,5) ,RYUM( 5,5) ,RUUM(5 ,5) ,RUYM(5 ,5) ,RUUINV(5,5) 
REAL  RYYINV(5,5),DM(5,5) ,CM( 5 ,5) .TRUNRM 

I NTEGER  RYYR , RYYC , RUUR , RUUC , RYUR , RYUC , RUYR , RUYC 
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DUA00170 

DUA00180 

DUA00190 

DUA00200 

DUA00210 

DUA00220 

DUA00230 

DUA00240 

DUA00250 

DUA00260 

DUA00270 

DUA00280 

DUA00290 

DUA00300 

DUA00310 

DUA00320 

DUA00330 

DUA00340 

DUA00350 

DUA00360 

DUA00370 

DUA00380 

DUA00390 

DUA00400 

DUA00410 

DUA00420 

DUA00430 

DUA00440 

DUA00450 

DUA00460 

DUA00470 

DUA00480 

DUA00490 

DUA00500 

DUA00510 

DUA00520 

DUA00530 

DUA00540 

DUA00550 

DUA00560 

DUA00570 

DUA00580 

DUA00590 

DUA00600 

DUA00610 

DUA00620 

DUA00630 

DUA00640 

DUA00650 

DUA00660 

DUA00670 
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INTEGER  DR,DC,CR,CC 

REAL  WKMAT( 12,12) ,X  TPO( 10, 10) ,Y( 2,2) ,U( 1) ,COEFM( 10, 10) 
INTEGER  WKR , WKC , XTR , XTC , COEFR , COEFC , ERR 

INTEGER  ITERA , CORRLN , YS I 2E , US I ZE , KTIME , SHFT 

DOUBLE  PRECISION  SEED 

*  TEMPORATY  MATRICES  FOR  PERFORMING  CALCUALTIONS 

REAL  T1M( 5,5), T2M( 5,5), T3M( 5,5) 

INTEGER  T1R , TIC , T2R , T2C , T3R , T3C 


*  BEGIN  MAIN  PROGRAM 


*  READ  IN  THE  SIZE  OF  THE  TEST  SYSTEM,  THE  NUMBER  OF  ITERATIONS 

*  TO  PERFORM  AND  THE  BEGINNING  LAG  OF  LINEAR  PREDICTION  OF  THE 

*  AR  PARAMETERS 

RE AD( 4 , * , END=22 )  YS I ZE , US I ZE , COEFR , COEFC , ITERA , SHFT 

*  READ  IN  THE  TRUE  VALUES  OF  THE  COEFFICIENTS  OF  THE  TEST  SYSTEM 

CALL  RDMATCCOEFM, COEFR .COEFC) 

*  INITIALIZE  ROW  AND  COLUMN  SIZES  OF  AS  WELL  AS  OTHER  VARIABLES 

CORRLN  =500 
RUUR  =  USIZE 
RUUC  =  USIZE 
RYYR  =  YSIZE  +  1 
RYYC  =  YSIZE  +  1 
T3R  =  YSIZE 
T3C  =  YSIZE 
DR  =  USIZE 
DC  =  1 

CR  =  YSIZE  +  1 
CC  =  1 
XTR  =  COEFC 
XTC  =  COEFR 
SEED  =  888. 0D1 
II  =  0 

Y(  1 , 1)  =  0.0 

E 1 ( 1,0)  =  0. 0 

DIVISR  =  TRUNRMCCOEFM, COEFR, YSIZE, USIZE) 

*  ZERO  OUT  THE  CORRELATION  MATRICES,  THE  MA  PARAMETER  VECTOR  AND 

*  VECTOR  OF  INPUT  DATA 

CALL  INITCRYYM, RYYR, RYYC, 0. 0) 

CALL  INIT( RUUM, RUUR, RUUC, 0.  0) 

CALL  INIT(RYUM, RYYR, RUUC, 0.  0) 

CALL  INIT( RUYM, RUUR, RYYC, 0.  0) 
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CALL  INIT(DM,DR,DC ,0. 0) 

CALL  INIT(X  TPO,XTR,XTC,0. 0) 

*  INITIALIZE  THE  AR  PARAMETER  VECTOR 

CALL  INITD(CM,CR,CC, 1.0) 

CALL  PRMAT( CM , CR , CC ) 

*  RUN  THE  FILTER  FOR  2000  TIME  STEPS  TO  GENERATE  OUTPUT  DATA 

jc - .... - - 


DO  70  KTIME  =  0,2000 

*  GET  VALUE  FOR  U(K).  U  IS  A  GAUSSIAN  RANDOM  VARIABLE. 

CALL  GGNML  (SEED,1,U) 

*  SHIFT  U(K)  INTO  X  VECTOR  TO  CREATE  X(K) 

CALL  SHIFT(X  TPO,XTR,XTC , YSIZE ,USIZE,Y( 1 , 1) ,U( 1) ) 

*  CALCULATE  VALUE  OF  Y(K) 

CALL  MULTI (X  TPO,XTR,XTC,COEFM,COEFR .COEFC ,Y, 1 , 1) 

*  SAVE  THE  INPUT  AND  OUTPUT  DATA 

RAWDATC1, KTIME)  =  Y(l,l) 

RAWDATC  2 , KTIME )  =  U(l) 

70  CONTINUE 

15  FORMAT  (2X,I4,2X,F8. 5,2X,F8.  5) 

*  CALCULATE  THE  CORRELATION  MATRICES 

CALL  AUTCOR ( 1 , RYYM , RYYR , RYYC , CORRLN ) 

CALL  AUTCOR( 2 , RUUM , RUUR , RUUC , CORRLN) 

CALL  CRSCOR( 1,2, RYUM , RYYR , RUUC , CORRLN ) 

CALL  CRSCOR( 2 , 1 , RUYM , RUUR , RYYC , CORRLN ) 

*  INVERT  THE  AUTOCORRELATION  MATRICES  OF  THE  OUTPUT  AND  INPUT  DATA 

CALL  LINV2F( RYYM , RYYR , RYYC ,RYYINV , 0 , WKMAT ,ERR) 

CALL  LINV2F( RUUM , RUUR , RUUC , RUUINV , 0 , WKMAT , ERR) 

*  MULTIPLY  THE  INVERSE  OF  THE  AUTOCORRELATION  MATRICES  BY  THEIR 

*  RESPECTIVE  CROSSCORRELATION  MATRICES 

CALL  MULTI (RUUINV, RUUR, RUUC, RUYM, RUUR, RYYC, TIM, T1R, TIC) 

CALL  MULTI (RYYINV, RYYR , RYYC , RYUM , RYYR , RUUC ,T2M ,T2R ,T2C ) 

*  ESTIMATE  THE  AR  PARAMETERS  BY  LINEAR  PREDICTION 

IF  (SHFT.GE.  1)  THEN 

CALL  C0RLA4( DR , CM , CR , CC ,T3M , CORRLN , SHFT) 

ENDIF 

WRITE  (3,26)  II ,DM( 1, 1) ,DM( 2,1) ,DM(3 , 1) ,DM(4, 1) ,DM(5 , 1) 

WRITE  (3,16)  II ,CM( 1,1) ,CM( 2,1) ,CM( 3, 1) ,CM(4, 1) ,CM(5 , 1) 

WRITE  (3,*)  ' 

CALL  SAVE( 1 , 0, II ,DM,DR,DC) 


DUA01240 

DUA01250 

DUA01260 

DUA01270 

DUA01280 

DUA01290 

DUA01300 

DUA01310 

DUA01320 

DUA01330 

DUA01340 

DUA01350 

DUA01360 

DUA01370 

DUA01380 

DUA01390 

DUA01400 

DUA01410 

DUA01420 

DUA01430 

DUA01440 

DUA01450 

DUA01460 

DUA01470 

DUA01480 

DUA01490 

DUA01500 

DUA01510 

DUA01520 

DUA01530 

DUA01540 

DUA01550 

DUA01560 

DUA01570 

DUA01580 

DUA01590 

DUA01660 

DUA01670 

DUA01680 

DUA01690 

DUA01780 

DUA01790 

DUA01800 

DUA01810 

DUA01820 

DUA01830 

DUA01840 

DUA01850 

DUA01870 

DUA01880 

DUA01890 

DUA01900 

DUA01910 

DUA01920 

DUA01930 

DUA01940 


1 


CALL  SAVEC  2 , 0 , I I , CM , CR , CC ) 


*  CALCULATE  THE  COEFFICIENT  ERROR 

CALL  ERR2 ( COEFM , COEFR , CM , CR , DM , DR , I I , YS I ZE , USI ZE , DI VI SR) 


*  BEGIN  THE  ITERATIVE  PROCEDURE 

*  _ _ ..... - - 


DO  100  II  =  1,ITERA 

*  CALCULATE  THE  VALUES  OF  THE  PARAMETERS  AT  ITERATION  II 

C ALL  MULTI ( T 1M , T1R , T 1C , CM , CR , CC , DM , DR , DC ) 

CALL  MULTI ( T2M , T2R , T2 C , DM , DR , DC , CM , CR , CC ) 

CALL  SAVE( 1,0, I I, DM, DR, DC) 

CALL  SAVE( 2 ,0, II ,CM ,CR,CC) 

*  CALCULATE  THE  STOPPING  PARAMETER 

CALL  ERR0R( CM , CR , CC , DM , DR , DC , I I ) 

*  CALCULATE  THE  COEFFICIENT  ERROR 

CALL  ERR2 ( COEFM , COEFR , CM , CR , DM , DR , I I , YS I ZE , US I ZE , D I VI SR ) 

WRITE  (3,26)  II ,DM( 1,1) ,DM(2 , 1) ,DM( 3, 1) ,DM(4, 1) ,DM(5 , 1) 

26  FORMAT  (I4,1X,'NUM  COEF  =’ ,5( IX, F9.  6)) 

WRITE  (3,1)  II ,CM( 1,1) ,CM( 2,1) ,CM( 3, 1) ,CM(4, 1) ,CM(5 , 1) 

16  FORMAT  (I4,1X,'DNM  COEF  =' ,5( 1X,F9. 6) ) 

WRITE  (3,*)  '  ERROR  =  ',£1(1,11),'  COEF  ERROR  =  ' ,E1(2,II) 
WRITE  (3,*)  ' 

CM(  1,1)  =  1.0 

100  CONTINUE 

*  END  OF  THE  ITERATIVE  PROCEDURE 


*  PLOT  THE  PARAMETER  VALUES 

CALL  PLOTKITERA, DR, CR, COEFM) 

22  STOP 
END 


* 


* 


FUNCTION  TRUNRM  (MAT1 ,M1R,YSZ ,USZ) 


REAL  MAT1(M1R, 1) 

INTEGER  YSZ.USZ 

VALUE  =  0 

DO  91  I  =  1 ,YSZ  +  USZ 

VALUE  =  VALUE  +  (MAT1( I , 1))**2 
91  CONTINUE 

TRUNRM  =  SQRT( VALUE  +  1) 

RETURN 
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END 


* 


* 


SUBROUTINE  AUTCOR ( IFST.MATl , MIR, MIC, CORRLN ) 


*  THIS  SUBROUTINE  CALCULATES  THE  AUTOCORRELATION  MATRIX  USING 

*  CORRELATIONS  OF  SIZE  CORRLN 


COMMON  /D /  RAWDAT(2,0: 2000) ,E1(2,0: 1000) , 

+ARRD(0:  1000,5) ,ARRN(0: 1000,5) 

REAL  MATI(MIR.MIC) 

INTEGER  I, J,K, CORRLN 

CALL  INIT(MAT1,M1R,M1C,0. 0) 

*  CALC  CORRELATIONS  ALONG  THE  FIRST  ROW  OF  THE  MATRIX.  THE 

*  SHIFT  =  ABS( NUMBER  OF  THE  COLUMN  -  1) 

*  LENGTH  OF  CORR  =  D  LENGTH  +  ABS(1  -  THE  NUMBER  OF  THE  COLUMN) 

*  ONCE  CORR  HAVE  BEEN  CALCULATED  FOR  THE  FIRST  ROW,  THEY  CAN 

*  BE  COPIED  INTO  OTHER  ROWS.  HORIZ  DISTANCE  OF  THE  PARTICULAR 

*  ELEMENT  FROM  THE  MAIN  DIAGONAL  DETERMINES  WHICH  CORR  TO 

*  COPY.  THIS  DISTANCE  IS  GIVEN  BY  ABS(ROW  NUMBER  -  COLUMN 

*  NUMBER). 

*  CORRELATIONS  START  200  POINTS  FROM  THE  BEGINNING  OF  THE 

*  DATA  TO  ELIMINATE  THE  TRANSIENT  OF  THE  TEST  SYSTEM. 

DO  90  J  »  1 ,M1C 

ENDPT  =  200  4-  CORRLN  -  ABS(1  -  J) 

DO  93  K  =  200, ENDPT 

MATl(l.J)  =  MAT1( 1 , J)  +  RAWDAT( IFST,K-( 1-J) )*RAWDAT( IFST 
93  CONTINUE 

MAT1(1,J)  =  MAT1( 1 , J)/(CORRLN  +  1  -  ABS(1  -  J)) 

90  CONTINUE 

DO  91  I  *  2, MIR 
DO  92  J  =  1 ,M1C 

MATl(I.J)  =  MAT1( 1 , 1+ABSC I-J) ) 

92  CONTINUE 

91  CONTINUE 

RETURN 

END 


* 


* 


SUBROUTINE  COPY( MAT 1 , M 1R , M 1C , MAT2 , M2R , M2C ) 


*  THIS  SUBROUTINE  COPIES  A  MATRIX  INTO  A  SECOND  MATRIX. 

REAL  MAT1(M1R,M1C),MAT2(M2R,M2C) 

INTEGER  I,J 

DO  90  I  ■  1 ,M1R 
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DO  900  J  =  1.M1C 

MAT2( I , J)  =  MAT1( I ,  J) 
900  CONTINUE 

90  CONTINUE 
M2R  =  MIR 
M2C  =  MIC 
RETURN 
END 


* 


* 


SUBROUTINE  CRSCORC IFST , ISND , MAT1 , MIR , MIC , CORRLN) 


*  THIS  SUBROUTINE  CALCULATES  THE  CROSSCORRELATION  MATRIX  USING 

*  CORRELATIONS  OF  LENGTH  CORRLN 


COMMON  /D/  RAWDAT( 2,0: 2000), El( 2,0: 1000), 

+ARRD( 0:  1000,5) ,ARRN( 0: 1000,5) 

REAL  MATI(MIR.MIC) 

INTEGER  I, J.K, CORRLN 

CALL  INIT(MAT1,M1R,M1C,0. 0) 

*  FOR  CROSS  CORRELATION  MUST  CALC  EACH  ELEMENT  OF  THE  CORR  MATRIX 

*  SEPARATELY.  CORRELATIONS  START  200  POINTS  FROM  THE  BEGINNING 

*  OF  THE  DATA  TO  ELIMINATE  THE  TRANSIENT  OF  THE  TEST  SYSTEM. 

DO  94  I  =  1 ,M1R 
DO  95  J  =  1 ,M1C 

ENDPT  =  200  +  CORRLN  -  ABS(I  -  J) 

IF  (J. GE. I)  THEN 

DO  96  K  =  200, ENDPT 

MAT1( I , J)  =  MATl(I.J)  +  RAWDAT( IFST,K-( I-J)) 
+*RAWDAT( ISND , K) 

96  CONTINUE 
ELSE 

DO  97  K  =  200, ENDPT 

MATl(I.J)  =  MATl(I.J)  +  RAWDAT(IFST.K)* 

+RAWDAT( 1 SND , K+( I - J ) ) 

97  CONTINUE 
END  IF 

MATl(I.J)  =  MAT1 ( I, J)/( CORRLN  +  1  -  ABS(I-J)) 

95  CONTINUE 

94  CONTINUE 
RETURN 
END 

SUBROUTINE  C0RLA4 ( IS I ZE , MAT2 , M2R , M2C , MAT3 , CORRLN , SHIFTT) 


*  THIS  SUBROUTINE  CALCULATES  THE  CORRELATION  MATRIX 

*  AND  THE  CORRELATION  VECTOR  USED  FOR  LINEAR  PREDICTION  OF  THE 

*  AR  PARAMETERS.  IT  THEN  CALCULATES  THE  INITIAL  ESTIMATE  OF  THE 

*  AR  PARAMETERS  AND  PASSES  THEM  BACK  TO  THE  MAIN  PROGRAM  IN  MAT2. 
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COMMON  /D/  RAWDAT( 2,0: 2000), El( 2,0: 1000), 

+ARRDC0:  1000,5) ,ARRN(0: 1000,5) 

REAL  MAT2(M2R,M2C) ,MAT3(M2R- 1 ,M2R-1) 

REAL  T2M( 5,1) ,T3M( 5,1) ,T4M( 5 ,5) ,WKMAT( 9,9) 

INTEGER  ERR.T1R ,T1C ,T2R ,T2C ,T3R,T3C , ISIZE ,SHIFTT,CORRLN 

*  GENERATE  THE  FIRST  ROW  OF  THE  RYY  MATRIX. 

*  SHIFTS  ARE  GREATER  THAN  THE  ORDER  OF  THE 

*  NUMERATOR. 

M3R  =  M2R-1 
M3C  =  M3R 

CALL  INIT(MAT3,M3R,M3C,0. 0) 

DO  90  J  =  1 ,M3C 

ENDPT  =  200  +  CORRLN  -  SHIFTT 
DO  91  K  =  200, ENDPT 

MAT3( 1 , J)  =  MAT3( 1 , J)  +  RAWDAT( 1,K+SHIFTT)*RAWDAT( 1,K) 

91  CONTINUE 

MAT3( 1 , J)  =  MAT3( 1 , J)/(C0RRLN-SHIFTT+1) 

SHIFTT  =  SHIFTT  +  1 

90  CONTINUE 

*  COPY  ELEMENTS  FROM  THE  FIRST  ROW  INTO  OTHER  LOCATIONS 

DO  92  I  =  2 ,M3R 
DO  93  J  =  l ,M3C 

MAT3( I , J)  =  MAT3( 1 , 1+ABS( I - J) ) 

93  CONTINUE 

92  CONTINUE 

*  GENERATE  THE  RYY  VECTOR  BY  COPYING  ELEMENTS  OF  THE  RYY  MATRIX 

T2R  =  M3C 
T2C  =  1 

CALL  FILL( 1 ,T2R- 1 , 1 , 1 ,T2M ,T2R , 1 , 2 , 1 ,MAT3 ,M3R ,M3C) 

*  GENERATE  THE  LAST  ELEMENT  IN  THE  RYY  CORRELATION  VECTOR. 
FINELE  =  0.0 

ENDPT  =  200  +  CORRLN  -  SHIFTT 
DO  94  K  =  200 .ENDPT 

FINELE  =  FINELE  +  RAWDATC 1 ,K+SHIFTT)*RAWDAT( 1,K) 

94  CONTINUE 

FINELE  *  FINELE/(C0RRLN+1 -SHIFTT) 

*  COPY  THE  FINAL  ELEMENT  INTO  THE  VECTOR 
T2M(T2R, 1)  =  FINELE 

*  CALCULATE  THE  INITIAL  ESTIMATE  OF  THE  AR  PARAMETERS 

CALL  LINV2F(MAT3,M3R,M3C,T4M,0,WKMAT,ERR) 

CALL  MULTI ( T4M , M3R , M3C , T2M , T2R , T2C , T3M , T3R , T3C ) 

*  COPY  THE  AR  PARAMETERS  INTO  THE  RETURN  ARGUMENT 
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( 


MAT2(  1,1)  =  1.0 
DO  95  L  =  1,3 

MAT2(L  +  1,1)  =  T3M( L, 1) 

*  MAT2(L, 1)  =  T3M(L- 1,1) 

*  WRITE  (3,*)  MAT2(L, 1) 

95  CONTINUE 

RETURN 

END 


SUBROUTINE  ERROR( MAT1, MIR, MIC, MAT2,M2R,M2C, ITNUM) 


*  THIS  SUBROUTINE  CALCULATES  THE  STOPPING  PARAMETER. 


COMMON  /D/  RAWDAT(2,0:  2000), El(2,0:  1000), 
+ARRD(0:  1000,5) ,ARRN(0: 1000,5) 

REAL  MAT1(M1R,M1C) ,MAT2(M2R,M2C) 

INTEGER  I , J,K 

ERVAL  =  0. 0 
XVAL  =  0. 0 
ZVAL  =  0.  0 

DO  90  I  =  400,450 
DO  91  J  =  MIR, 1,-1 

XVAL  =  XVAL  +  MAT1( J, 1)*RAWDAT( 1 , I+M1R-J) 

91  CONTINUE 

DO  92  K  =  M2R ,1,-1 

ZVAL  =  ZVAL  +  MAT2(K, 1)*RAWDAT( 2 , I+M2R-K) 

92  CONTINUE 

ERVAL  =  ERVAL  +  (XVAL-ZVAL)**2 
XVAL  =  0.  0 
ZVAL  =  0.  0 
90  CONTINUE 

E 1 ( 1, ITNUM)  =  ERVAL/51 

RETURN 

END 


* 


* 


Meiti 


SUBROUTINE  ERR2 ( MAT 1 , MIR , MAT 2 , M2R , MAT 3 , M3R , ITNUM , 
+YSZ ,USZ,DIVISR) 


*  THIS  SUBROUTINE  CALCUALTES  THE  COEFFICIENT  ERROR. 

COMMON  /D/  RAWDAT(2,0: 2000), El(2,0: 1000), 
+ARRD(0:  1000,5) ,ARRN(0: 1000,5) 

REAL  MATKM1R,  1)  ,MAT2(M2R,  1)  ,MAT3(M3R,  1) 
INTEGER  I,YSZ,USZ 

ERRVAL  =  (1  -  MAT2(1,1))**2 
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DO  90  I  =  l.YSZ 

ERRVAL  =  ERRVAL  +  (MAT1(I,1)  +  MAT2( 1+1 , 1) )**2 
90  CONTINUE 

DO  92  I  =  l.USZ 

ERRVAL  =  ERRVAL  +  (MAT1(I+YSZ,1)  -  MAT3( I , 1) )**2 
92  CONTINUE 


E1(2,ITNUM)  =  SQRT( ERRVAL) /DIVISR 

RETURN 

END 


* 


* 


SUBROUTINE  FILL(I,J,K,L,MAT1,M1R,M1C,R2,C2,MAT2,M2R,M2C) 


*  THIS  ROUTINE  FILLS  MAT1  FROM  MAT2.  POSITIONS 

*  IN  MAT1  FROM  (I,K)  TO  (J,L)  ARE  FILLED  WITH  AN  EQUAL  NUMBER 

*  OF  ELEMENTS  FROM  MAT2  STARTING  AT  POSITION  (R2,C2). 

REAL  MAT1(M1R,M1C) ,MAT2(M2R,M2C) 

INTEGER  ROW,COL,R2,C2 ,C22 

C22  =  C2 

DO  90  ROW  =  I,J 
DO  91  COL  =  K,L 

MATl(ROW,COL)  =  MAT2(R2,C2) 

C2  =  C2  +  1 
91  CONTINUE 

R2  =  R2  +  1 
C2  —  C22 
90  CONTINUE 


RETURN 

END 

SUBROUTINE  LI M2 ( EMAX, ESTP.CEMAX.CESTP, ITERA) 


*  ROUTINE  CALCULATES  THE  LIMITS  FOR  THE  GRAPH  OF  THE  STOPPING 

*  PARAMETER  AND  THE  COEFFICIENT  ERROR. 

*  STOPPING  PARAMETER  LIMITS  ARE  RETURNED  IN  EMAX  AND  ESTP. 

*  COEFFICIENT  ERROR  LIMITS  ARE  RETURNED  IN  CEMAX  AND  CESTP. 

COMMON  /D/  RAWDAT(2,0: 2000), El(2,0: 1000), 

+ARRDC0:  1000,5) ,ARRN(0:  1000,5) 

REAL  EMAX, ESTP 
INTEGER  ITERA 


EMAX  =  0.  0 

DO  90  I  =  0, ITERA 

IF  (El( 1,1). GT. EMAX)  THEN 
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EMAX  =  E 1 ( 1,1)  DUA07190 

ENDIF  DUA07200 

90  CONTINUE  DUA07220 

EMAX  =  1.  25  *  EMAX  DUA07230 

ESTP  =  EMAX/5  DUA07240 

CEMAX  =0.0  DUA07280 


DO  91  I  =  O.ITERA  DUA07290 

IF  (E1(2,I).GT. CEMAX)  THEN  DUA07310 

CEMAX  =  E1(2,I)  DUA07320 

ENDIF  DUA07330 

91  CONTINUE  DUA07350 


CEMAX  =  1.  25  *  CEMAX 
CESTP  =  CEMAX/ 5 
RETURN 
END 

SUBROUTINE  PLOT 1 ( ITERA , I CUR VN , I CURVD , MAT 1 ) 
*  ****************************** 


*  THIS  ROUTINE  GENERATES  SEPARATE  PLOTS  OF  THE  MA 

*  AND  AR  PARAMETERS.  IT  THEN  REPLOTS  THE  THESE  CURVES  ALONG 

*  WITH  THE  STOPPING  PARAMETER.  FINALLY  IT  PLOTS  THE  STOPPING 

*  PARAMETER  AND  THE  COEFFICIENT  ERROR  ON  THE  SAME  GRAPH. 

COMMON  /D/  RAWDAT(2,0:  2000), £1(2,0: 1000), 

+ARRD( 0:  1000 ,5) ,ARRN( 0:  1000,5) 

REAL  X(0: 1000), Y(0: 1000) ,MAT1( 10, 10) 

REAL  STP,RNMIN,RNSTEP,RNMAX,RITERA 
INTEGER  I , J , I 1R , I 1C , VAL , ITERA , ICURV 

CALL  LIMITS ( ICURVN, ICURVD,DMAX,DMIN,DSTEP, 

+NMAX , NMIN , NSTEP , ITERA ) 

CALL  LIM2( EMAX , ESTP , CEMAX , CESTP , ITERA ) 

*  GENERATE  THE  ITERATION  NUMBER 

DO  90  I  =  O.ITERA 
X(I)  =  I 
Y( I)  =  0.  0 
90  CONTINUE 

*  CALCULATE  X  AX  S  JELING  INTERVAL 

STP  =  ITER.  ,  10. 

*  SET  THE  DEVICE  CALLS  FOR  ALL  OF  THE  PLOTS 

CALL  SHERPA( 'MCGRAF  ','A',3) 

CALL  RESET(’ALL') 

*  SECTION  1 

*  THIS  SECTION  GRAPHS  THE  NUMERATOR  AND  DENOMINATOR 

*  PARAMETERS  ON  SEPARATE  GRAPHS  ON  THE  SAME  PAGE. 
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DUA07930 

DUA07940 

DUA07950 

DUA07970 

DUA07980 

DUA07990 

DUA08000 

DUA08030 

DUA08050 

DUA08060 

DUA08070 

DUA08080 

DUA08090 
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*  SET  UP  THE  PLOT  FOR  THE  NUMERATOR  COEFF 

CALL  PAGE( 8. 5,11. 0) 

CALL  HEIGHTCO.  14) 

CALL  HWROT( ’AUTO* ) 

CALL  XINTAX 

CALL  PHYS0R(1. 5,6.  0) 

CALL  AREA2D(5. 0,3. 5) 

CALL  COMPLX 
CALL  YAXANG(O) 

CALL  XNAMEC ’ ITERATIONS? ' , 100) 

CALL  YNAMEC 'PARAMETER  VALUES? ', 100) 

CALL  MESSAG('(A)$' ,100,2.  4,-0. 8) 

CALL  THKFRM(0. 03) 

CALL  FRAME 

CALL  GRAF(0.  , STP , ITERA , NMIN , NSTEP , NMAX) 

*  PLOT  THE  NUMERATOR  PARAMETERS 

DO  93  J  =  1 , ICURVN 
DO  94  I  =  0, ITERA 

Y( I)  =  ARRN(I,J) 

94  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

93  CONTINUE 

*  PLOT  DASHED  LINES  FOR  THE  TRUE  VALUE  OF  THE  PARAMETERS 

CALL  DASH 

DO  97  K  =  ICURVD, ICURVD+ICURVN-1 
DO  98  J  *  0, ITERA 
Y(J)  =  MAT1(K, 1) 

98  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

97  CONTINUE 

CALL  ENDGR(O) 

*  SET  UP  THE  PLOT  FOR  THE  DENOMINATOR  PARAMETERS 

CALL  RESETC ’DASH’ ) 

CALL  PHYSOR( 1.5,1.  5) 

CALL  AREA2D(5.  0,3.  5) 

CALL  XNAMEC ' ITERATIONS? ' ,100) 

CALL  YNAMEC ' PARAMETER  VALUES? 100) 

CALL  MESSAGC '(B)?' , 100,2.  4,-0.  8) 

CALL  THKFRMCO.  03) 

CALL  FRAME 

CALL  GRAF(0. , STP, ITERA, DMIN,DSTEP,DMAX) 

*  PLOT  THE  DENOMINATOR  PARAMETERS 

DO  95  J  =  1, ICURVD 
DO  96  I  =  0, ITERA 

Y( I)  =  ARRDC I , J) 

96  CONTINUE 


DUA08100 

DUA08110 

DUA08120 

DUA08130 

DUA08140 

DUA08150 

DUA08160 

DUA08170 

DUA08200 

DUA08210 

DUA08240 

DUA08250 

DUA08260 

DUA08270 

DUA08280 

DUA08290 

DUA08320 

DUA08340 
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DUA08410 

DUA08420 

DUA08430 

DUA08440 

DUA08450 

DUA08460 

DUA08470 

DUA08480 

DUA08490 

DUA08500 

DUA08510 

DUA08520 

DUA08530 

DUA08540 

DUA08550 

DUA08560 

DUA08570 

DUA08590 

DUA08600 

DUA08620 

DUA08630 

DUA08640 

DUA08650 

DUA08660 

DUA08680 

DUA08700 

DUA08710 

DUA08720 

DUA08730 

DUA08740 

DUA08750 

DUA08760 
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I 


CALL  CURVE(X,Y, ITERA, 0) 

95  CONTINUE 

*  PLOT  DASHED  LINES  FOR  THE  TRUE  VALUES  OF  THE  DENOM  PARAMETERS 

CALL  DASH 

DO  99  K  »  1 , ICURVD-1 
DO  100  J  =  O.ITERA 
Y(J)  =  -MATl(K.l) 

100  CONTINUE 

CALL  CURVE( X , Y , ITERA , 0 ) 

99  CONTINUE 

DO  105  J  =  O.ITERA 
Y(J)  *1.0 
105  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

CALL  ENDPL(O) 

*  SECTION  2 

*  THIS  SECTION  PUTS  THE  STOPPING  PARAMETER  ON  THE 

*  GRAPHS  OF  THE  NUMERATOR  AND  DENOMINATOR  PARAMETERS. 

*  SET  UP  THE  PLOT  FOR  THE  NUMERATOR  PARAMETERS 

CALL  RESET('DASH') 

CALL  HWROT('AUTO') 

CALL  XINTAX 

CALL  PHYSOR(l.  5,6.0) 

CALL  AREA2DC5.  0,3.  5) 

CALL  COMPLX 
CALL  YAXANG(O) 

CALL  XNAME( ' ITERATIONS$ ' , 100) 

CALL  YNAMEC' PARAMETER  VALUES*' , 100) 

CALL  MESSAG( ' ( A) $ ’ ,100,2.  4,-0.  8) 

CALL  THKFRM(0.  03) 

CALL  FRAME 

CALL  GRAF(0. ,STP, ITERA, NMIN,NSTEP,NMAX) 

*  PLOT  THE  NUMERATOR  PARAMETERS 

DO  200  J  ~  1 , ICURVN 
DO  201  i  *  O.ITERA 

Y( I)  =  ARRN(I,J) 

201  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

200  CONTINUE 

*  PLOT  DASHED  LINES  FOR  THE  TRUE  VALUES  OF  THE  PARAMETERS 

CALL  DASH 

DO  202  K  =  ICURVD , ICURVD+ICURVN- 1 
DO  203  J  =  O.ITERA 
Y( J)  =  MAT1(K,1) 

203  CONTINUE 


DUA08770 

DUA08780 

DUA08790 

DUA08800 

DUA08810 

DUA08820 

DUA08830 

DUA08840 

DUA08850 

DUA08860 

DUA08870 

DUA08880 

DUA08890 

DUA08900 

DUA08910 

DUA08920 

DUA08930 

DUA08940 

DUA08950 

DUA08960 

DUA08970 

DUA08980 

DUA08990 

DUA09010 

DUA09020 

DUA09030 

DUA09040 

DUA09050 

DUA09060 

DUA09090 

DUA09100 

DUA09130 

DUA09140 

DUA09150 

DUA09160 

DUA09170 

DUA09180 

DUA09190 

DUA09210 

DUA09220 

DUA09230 

DUA09240 

DUA09250 

DUA09260 

DUA09270 

DUA09280 

DUA09290 

DUA09300 

DUA09310 

DUA09320 

DUA09330 

DUA09340 

DUA09350 

DUA09360 

DUA09370 
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CALL  CURVE (X,Y, ITERA, 0) 
202  CONTINUE 


*  PLOT  THE  STOPPING  PARAMETER  ON  THE  SAME  GRAPH 
CALL  DOT 

CALL  YGRAXS(0.  0,ESTP,EMAX,3. 5, 'STOPPING  PARAMETER?', 
+-100,5. 0,0. 0) 

DO  204  J  =  0, ITERA 
Y(J)  =  El( 1 , J) 

204  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

CALL  ENDGR(O) 


*  SET  UP  THE  PLOT  FOR  THE  DENOMINATOR  PARAMETERS 

CALL  RF.SET('DOT') 

CALL  PHYSOR( 1. 5 , 1. 5) 

CALL  AREA2D(5. 0,3.5) 

CALL  XNAME( ' ITERATIONS? ' , 100) 

CALL  YNAME(' PARAMETER  VALUES? ’, 100) 

CALL  MESSAG( ' (B)? 1 , 100,2.  4, -0. 8) 

CALL  THKFRM( 0. 03) 

CALL  FRAME 

CALL  GRAF(0.  ,STP, ITERA, DMIN,DSTEP,DMAX) 

*  PLOT  THE  DENOMINATOR  PARAMETERS 

DO  205  J  =  1 , ICURVD 
DO  206  I  =  0, ITERA 

Y( I )  =  ARRD( I , J) 

206  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

205  CONTINUE 

*  PLOT  DASHED  LINES  FOR  THE  TRUE  VALUES  OF  THE  PARAMETERS 

CALL  DASH 

DO  207  K  =  1, ICURVD- 1 
DO  208  J  =  0, ITERA 
Y( J)  =  -MAT1(K, 1) 

208  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

207  CONTINUE 

DO  209  J  =  0, ITERA 
Y(  J)  =  1.0 

209  CONTINUE 

CALL  CURVE(X,Y, ITERA, 0) 

*  PLOT  THE  STOPPING  PARAMETER  ON  THE  SAME  GRAPH 

CALL  DOT 


DUA09380 

DUA09390 

DUA09400 

DUA09410 

DUA09420 

DUA09430 

DUA09440 

DUA09450 

DUA09460 

DUA09470 

DUA09480 

DUA09490 

DUA09500 

DUA09510 

DUA09520 

DUA09530 

DUA09540 

DUA09550 

DUA09560 

DUA09570 

DUA09580 

DUA09590 

DUA09610 

DUA09620 

DUA09630 

DUA09640 

DUA09650 

DUA09660 

DUA09680 

DUA09690 

DUA09700 

DUA09710 

DUA09720 

DUA09730 

DUA09740 

DUA09750 

DUA09760 

DUA09770 

DUA09780 

DUA09790 

DUA09800 

DUA09810 

DUA09820 

DUA09830 

DUA09840 

DUA09850 

DUA09860 

DUA09870 

DUA09880 

DUA09890 

DUA09900 

DUA09910 

DUA09920 

DUA09930 

DUA09940 

DUA09950 


CALL  YGRAXSCO.  0,ESTP,EMAX, 3. 5, 'STOPPING  PARAMETER? ' , 
+-100,5. 0,0. 0) 


DO  210  J  =  0 , 1TERA 
Y(J)  =  El( 1 , J) 

210  CONTINUE 

CALL  CURVE (X, Y, ITERA,0) 

CALL  ENDPL(O) 

*  SECTION  3 

*  THIS  SECTION  PLOTS  THE  STOPPING  PARAMETER  AND  THE  COEFFICIENT 

*  ERROR  ON  THE  SAME  GRAPH. 

*  SETUP  THE  PLOT  FOR  THE  STOPPING  PARAMETER 

CALL  DOT 

CALL  HWROT('AUTO’) 

CALL  PHYSORU.  5,6.0) 

CALL  AREA2D(5. 0,3.  5) 

CALL  XNAME( ' ITERATIONS? ' , 100) 

CALL  YNAME( ' STOPPING  PARAMETER? ', 100) 

CALL  THKFRM(0.  03) 

CALL  FRAME 

CALL  GRAF( 0.  ,STP, ITERA.O.  ,ESTP,EMAX) 


DO  306  J  =  0 , ITERA 
Y(J)  =  El( 1 , J) 

306  CONTINUE 

CALL  CURVE(X,Y, ITERA.O) 


*  PLOT  THE  COEFFICIENT  ERROR  ON  THE  SAME  GRAPH 
CALL  CHNDOT 

CALL  YGRAXSCO.  0 ,CESTP ,CEMAX, 3. 5,’ COEFFICIENT  ERROR?’, 
+-100,5. 0,0. 0) 


DO  307  J  =  0, ITERA 
Y( J)  =  El( 2 , J) 

307  CONTINUE 

CALL  CURVE(X,Y, ITERA.O) 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 

SUBROUTINE  SAVE(VAL,M,K,MAT1,M1R,M1C) 


*  ROUTINE  SAVES  PARAMETER  ESTIMATES  IN  EITHER  ARRD  OR  ARRN 

*  DEPENDING  UPON  THE  VALUE  OF  VAL. 

COMMON  /D/  RAWDAT(2,0: 2000), El( 2,0: 1000), 

+ARRDCO: 1000,5) ,ARRN(0: 1000,5) 


DUA09960 
DUA09970 
DUA09980 
DUA09990 
DUA10000 
DUA10010 
DUA 10020 
DUA 100 30 
DUA 10 040 
DUA10050 
DUA10060 
DUA10070 
DUA 10080 
DUA10090 
DUA10100 
DUA10110 
DUA10120 
DUA10130 
DUA10140 
DUA10150 
DUA10160 
DUA10180 
DUA10190 
DUA10200 
DUA10220 
DUA10230 
DUA10240 
DUA10250 
DUA10260 
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DUA10280 
DUA10290 
DUA10300 
DUA10310 
DUA10320 
DUA10330 
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DUA10360 
DUA10370 
DUA10380 
DUA10390 
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DUA10890 
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DUA10940 


67 


REAL  MAT1(M1R,M1C) 

DUA10960 

INTEGER  I,J,K,M,VAL 

DUA10970 

DUA 10980 

DO  90  I  =  1,M1R 

DUA10990 

DO  900  J  =  1,M1C 

DUA 11000 

IF  (VAL. EQ. 1)  THEN 

DUA11010 

ARRN(K, I)  =  MAT1( I , J) 

DUA11020 

ELSE IF  (VAL.  EQ.  2)  THEN 

DUA11030 

ARRD(K,I)  =  MATl(I.J) 

DUA11040 

ENDIF 

DUA11050 

CONTINUE 

DUA11060 

CONTINUE 

DUA11070 

DUA11080 

RETURN 

DUA11090 

END 

DUA11100 

C.  SUBPROGRAMS  COMMON  TO  BOTH  SYSTEM  IDENTIFICATION 
ALGORITHMS  SUB00010 


* 

* 


920 

92 


SUB00020 

SUB00040 
SUB00050 
SUB00060 
SUB00070 

THIS  SUBROUTINE  ADDS  TWO  EQUAL  SIZE  MATRICIES  AND  PUTS  THE  RESULT  SUB00080 


SUBROUTINE  ADD  (MAT1 , IR1 , IC1 ,MAT2 , IR2 , IC2 ,RMAT, IRR, IRC) 


95 

94 


IN  A  THIRD  MATRIX. 

REAL  MAT1( IR1 , IC1) ,MAT2( IR2 , IC2) ,RMAT( IR1, IC1) 
INTEGER  I, J, IRR, IRC 

DO  92  1=1, IR1 
DO  920  J=1,IC1 

RMAT( I , J)  =  MAT1( I , J)  +  MAT2(I,J) 
CONTINUE 
CONTINUE 
IRR  =  IR1 
IRC  =  IC1 
RETURN 
END 

SUBROUTINE  1NIT( MAT1 ,M1R ,M1C , INITVL) 


THIS  SUBOUTINE  INITIALIZES  A  MATRIX  TO  INITVL 

REAL  MATKM1R, MIC), INITVL 
INTEGER  I,J 

DO  94  1=1, MIR 
DO  95  J=1,M1C 

MAT1( I , J)=INITVL 
CONTINUE 
CONTINUE 
RETURN 


SUB00090 

SUB00100 

SUB00130 

SUB00140 

SUB00160 

SUB00170 

SUB00180 

SUB00190 

SUB00200 

SUB00210 

SUB00220 

SUB00230 

SUB00240 

SUB00250 

SUB00260 

SUB00280 

SUB00290 

SUB00300 

SUB00310 

SUB00320 

SUB00330 

SUB00340 

SUB00350 

SUB00360 

SUB00380 

SUB00390 

SUB00400 

SUB00410 

SUB00420 

SUB00430 


* 

* 

* 

* 


95 

94 


* 

* 

* 

* 

* 


* 


91 

90 


END 

★★■ft*************************** ****** 

SUBROUTINE  INITD(MAT1 , MIR, MIC, INITVL) 
************************************ 


THIS  SUBROUTINE  INITIALIZES  A  MATRIX  TO  AS  A  DIAGONAL  MATRIX 
WHOSE  DIAGONAL  ELEMENTS  EQUAL  INITVL. 

REAL  MATKM1R, MIC), INITVL 
INTEGER  I,J 

DO  94  1=1, MIR 
DO  95  J=1 ,M1C 

IF  (I.EQ.  J)  THEN 
MAT1( I , J)=INITVL 
ELSE 

MAT1(  I ,  J)=0.  0 
END  IF 
CONTINUE 
CONTINUE 
RETURN 
END 


SUBROUTINE  LIMITS(NSZ,DSZ,DMAXjDMIN,DSTEP,NMAX,NMIN,NSTEP,ITERA) 


ROUTINE  CALCULATES  THE  LIMITS  FOR  THE  GRAPHS 
CALCULATES  DENOMINATOR  AND  NUMERATOR  LIMITS  SEPARATELY 
IN  PREPARATION  FOR  MAKING  TWO  GRAPHS 

COMMON  /D/  RAWDAT(2,0: 2000), El(2,0: 1000), 

+ARRD( 0: 1000,5) ,ARRN(0: 1000,5) 

REAL  DMAX,DMIN,DSTEP,NMAX,NMIN,NSTEP 
INTEGER  DSZ.NSZ 

CALCULATE  THE  DENOMINATOR  LIMITS 


DMAX  =1.0 
DMIN  =  0.0 

DO  90  I  =  1 ,DSZ 

DO  91  J  =  0 , ITERA 
IF  ( (ARRD( J, I)).  GT.  DMAX)  THEN 
DMAX  =  ARRD( J, I) 

ENDIF 

IF  ((ARRD(J,I)).  LT.DMIN)  THEN 
DMIN  =  ARRD( J, I) 

ENDIF 

CONTINUE 

CONTINUE 

IF  (DMAX. GT. 0)  THEN 
DMAX  =  1.  25  *  DMAX 
ELSE 


SUB00440 
SUB00450 
SUB00470 
SUB00480 
SUB00490 
SUB00500 
SUB005 10 
SUB00520 
SUB00530 
SUB00540 
SUB00550 
SUB00570 
SUB00580 
SUB00590 
SUB00600 
SUB00610 
SUB00620 
SUB00630 
SUB00640 
SUB00650 
SUB00660 
SUB00670 
SUB00680 
SUB00690 
SUB00710 
SUB00720 
SUB00730 
SUB00740 
SUB00750 
SUB00760 
SUB00770 
SUB00780 
SUB00790 
SUB00800 
SUB00820 
SUB00830 
SUB00840 
SUB00850 
SUB00860 
SUB00870 
SUB00880 
SUB00890 
SUB00900 
SUB00910 
SUB00920 
SUB00930 
SUB00940 
SUB00950 
SUB00960 
SUB00970 
SUB00980 
SUB00990 
SUB01000 
SUB01010 
SUB01020 
SUB01030 
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DMAX  =  0.0 
END  IF 

IF  (DMIN.GT.  0)  THEN 
DMIN  =  0.0 
ELSE 

DMIN  =  1.  25  *  DMIN 
ENDIF 

DSTEP  =  (DMAX  -  DMIN)/5 

CALCULATE  THE  NUMERATOR  LIMITS 

NMAX  =  0.  0 
NMIN  =  0.  0 
DO  92  I  =  l.NSZ 

DO  93  J  =  O.ITERA 
IF  (ARRN(J.I).GT.NMAX)  THEN 
NMAX  =  ARRN( J, I ) 

ENDIF 

IF  (ARRN(J.I).LT.  NMIN)  THEN 
NMIN  =  ARRN( J, I) 

ENDIF 

CONTINUE 

CONTINUE 

IF  (NMAX. GT. 0)  THEN 
NMAX  =  1.  25  *  NMAX 
ELSE 

NMAX  =  0.  0 
ENDIF 

IF  (NMIN.  GT.  0)  THEN 
NMIN  =  0.0 
ELSE 

NMIN  =  1. 25  *  NMIN 
ENDIF 

NSTEP  =  ABS(NMAX  -  NMIN)/5 

RETURN 

END 

ft&iv&lc'kit&irit&irftftjrirk'fr'trlc'ir'k'k'k'k'frkir'k'k'k'/c'tc'k’tclt'k'b’ft'k'k'k'tck'k'k'k+nk'k 

SUBROUTINE  MULTI  (MAT1,M1R,M1C,MAT2,M2R,M2C,RMAT,M3R,M3C) 


ROUTINE  MULTIPLIES  TWO  MATRICES  AND  PUT  THE  RESULT  IN  A 
THIRD  MATRIX. 

REAL  MATI(MIR.MIC) ,MAT2(M2R,M2C) ,RMAT(M1R,M2C) 

INTEGER  I , J,K, IRR, IRC 


CALL  INIT(RMAT,M1R,M2C,0.  0) 

DO  91  1=1, MIR 
DO  910  J=1 ,M2C 


SUB01040 
SUB01050 
SUB01060 
SUB01070 
SUBO 1080 
SUB01090 
SUB01100 
SUB01110 
SUB01130 
SUB01140 
SUB01150 
SUB01160 
SUB01170 
SUB01180 
SUB01190 
SUB01200 
SUB01210 
SUB01220 
SUBO 1230 
SUB01240 
SUB01250 
SUB01260 
SUB01270 
SUB01280 
SUB01290 
SUB01300 
SUB01310 
SUB01320 
SUB01330 
SUB01340 
SUB01350 
SUB01360 
SUB01370 
SUB01380 
SUB01390 
SUB01400 
SUB01410 
SUBO 1430 
SUB01440 
SUB01450 
SUBO 1460 
SUB01470 
SUB01490 
SUB01500 
SUB01510 
SUB01520 
SUB01530 
SUB01540 
SUBO 1550 
SUB01580 
SUB01590 
SUB01600 
SUB01610 
SUB01620 
SUB01630 
SUB01640 
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DO  9100  K=1,M1C 

RMAT( I , J)=RMAT( I , J)  +MAT1( I ,K)*MAT2(K, J) 
9100  CONTINUE 

910  CONTINUE 

91  CONTINUE 
M3R  =  MIR 
M3C  =  M2C 
RETURN 
END 


SUBROUTINE  PRMAT(MAT1,I1R,I1C) 


SUB01650 

SUB01660 

SUB01670 

SUB01680 

SUB01690 

SUB01700 

SUB01710 

SUB01720 

SUB01730 

SUB01740 

IVA04910 

IVA04920 

IVA04930 


*  SUBROUTINE  PRINTS  A  MATRIX  OUT  TO  THE  FILE  DEFINED 

*  AS  UNIT  3 

I  REAL  MAT1( 10, 10) 

INTEGER  I,J 

DO  92  I  =  1,I1R 

WRITE( 3,302)  (MAT1(I,J),J  =  1,I1C) 

302  FORMAT  (7(2X,F8.5)) 

j  92  CONTINUE 

RETURN 

END 

*  ^*************************** 

SUBROUTINE  RDMAT( MAT1 ,M1R,M1C) 

I  , 

*  ROUTINE  READS  A  MATRIX  FROM  FILE  SPECIFIED  AS  UNIT  4. 

REAL  MAT1(M1R,M1C) 

INTEGER  I,J 

*  READ  IN  MATRIX 

DO  92  I  *  1,M1R 

READ  (4,*)  (MAT1( I , J ) ,  J=1 ,M1C) 

C301  FORMATC 10F3. 1) 

92  CONTINUE 

RETURN 
END 


SUBROUTINE  SHIFTCMATI.MIR.MIC.DSIZE.NSIZE.OUTDAT.INDAT) 


*  ROUTINE  SHIFTS  NEW  INPUT  AND  OUTPUT  VALUES  INTO  DATA  VECTOR 

*  OF  THE  TEST  SYSTEM.  THE  OLDEST  VALUES  ARE  LOST  TO  MAKE 

*  ROOM  FOR  THE  NEW  VALUES. 

REAL  MAT1(M1R,M1C) ,OUTDAT( 1) ,INDAT( 1) 

INTEGER  J,DSIZE ,NSIZE ,NSTART 


IVA04970 

1VA04980 

IVA04990 

IVA05010 

IVA05020 

IVA05030 

IVA05040 

IVA05050 

IVA05060 

IVA05070 

SUB01760 

SUB01770 

SUB01780 

SUB01790 

SUB01800 

SUB01810 

SUB01840 

SUB01850 

SUB01860 

SUB01870 

SUB01880 

SUB01890 

SUB01900 

SUB01910 

SUB01920 

SUB01930 

SUB01940 

SUB01950 

SUB01970 

SUB01980 

SUB01990 

SUB02000 

SUB02010 

SUB02020 

SUB02030 

SUB02060 

SUB02070 

SUB02080 

SUB02090 
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NSTART  =  DSIZE  +  1 

SUB02100 

DO  92  J  =  DSIZE, 2,-1 

SUB02110 

MAT1(1,J)  =  MAT1( 1, J-l) 

SUB02120 

92 

CONTINUE 

SUB02130 

MAT1( 1,1)  =  OUTDAT(l) 

SUB02140 

SUB02150 

DO  93  J  =  MIC ,NSTART+1 , - 1 

SUB02160 

MAT1( 1 , J)  =  MAT1( 1 , J-l) 

SUB02170 

93 

CONTINUE 

SUB02180 

MAT 1(1, NSTART)  =  INDAT(l) 

SUB02190 

WRITE( 12,19)  (MAT1( 1 , J) , J=1 ,M1C) 

SUB02200 

19 

FORMAT  ( 7( IX, F8.  4) ) 

SUB02210 

SUB02220 

RETURN 

SUB02230 

END 

SUB02240 

SUB02260 

* 

irk*************K*irk**-;rt<-ki<-ir****ir-k*lrkirk**-irir* 

SUB02270 

SUBROUTINE  SMULTI  (CONST, MAT1 , MIR, MIC) 

SUB02280 

* 

SUB02290 

SUB02300 

* 

ROUTINE  MULTIPLIES  A  MATRIX  BY  A  SCALAR. 

SUB02310 

SUB02320 

REAL  MAT1(M1R, MIC), CONST 

SUB02350 

INTEGER  I,J 

SUB02370 

SUB02380 

DO  93  1=1, MIR 

SUB02390 

DO  930  J=1 ,M1C 

SUB02400 

MAT1( I , J)  =  MAT1( I , J)  *  CONST 

SUB02410 

930 

CONTINUE 

SUB02420 

93 

CONTINUE 

SUB02430 

SUB02440 

RETURN 

SUB02450 

END 

SUB02460 

SUB02480 
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